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slender elliptic cone as a model for non-linear 
supersonic flow theory 


By M. D. VAN DYKE 


Ames Aeronautical Laboratory, National Advisory Committee for Aeronautics, 
Moffett Field, California 


(Received 12 December 1955) 


SUMMARY 
The second-order slender-body solution is derived for an 
unyawed elliptic cone in supersonic flow. ‘The result is used as 
the basis for a critique of various approximations in compressible 
flow theory: slender-body, linearized, first- and second-order 
thin-wing theories ; edge corrections ; and the method of linearized 
characteristics. 


1. INTRODUCTION 

‘The circular cone serves as a standard of comparison for supersonic 
flow past bodies of revolution, and the elliptic cone will probably assume 
the same role for non-circular shapes. ‘Exact’ numerical solutions for 
elliptic cones comparable with those available for circular cones would 
require elaborate computing programmes that are not yet contemplated. 
However, a number of approximate theories for the elliptic cone have 
recently been advanced, based upon various linearizing assumptions. It 
can be anticipated that attempts will be made to improve some of these by 
successive approximations, so as to take account of the non-linear nature 
of compressible flow. 

In this paper a second approximation to the supersonic flow past un- 
yawed elliptic cones is obtained by proceeding from the slender-body 
theory. ‘The first approximation, following Ward’s (1949) definitive 
treatment of supersonic slender-body theory, has been given by Fraenkel 
(1952). ‘The principles of deriving the second approximation by iteration 
have been set forth by Lighthill (1954) and Van Dyke (1956), the present 
treatment of elliptic cones being the first application to non-circular bodies. 
Whereas Adams & Sears’ ‘ not-so-slender-body ’ theory (1953) seeks only 
a closer approach to the solution of the linearized Prandtl-Glauert equation, 
second-order slender-body theory includes also the leading non-linear terms. 

Although the solution given here has intrinsic interest, it is utilized 
primarily as a model for the full inviscid solution, and thus serves as the 
basis for a critique of various approximations commonly employed ip 
compressible flow theory. ‘Thus the second-order solution is regarded as 
being exact (since in one respect or another it is indeed more exact than any 
of the approximations to be considered), and the approximation to be tested 
is introduced in addition to the simplifications already made. 


F.M. 
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2. THE SECOND-ORDER SLENDER-BODY SOLUTION 
2.1. Résumé of the first-order solution 
Consider a uniform supersonic stream of Mach number / flowing along 
the z-axis. Introduce the non-orthogonal elliptic-conical coordinates 
(é, s) by setting 
x+1y =cs cosh (€ +7), 


(1) 


Then an unyawed elliptic cone is described by é = &,j or (figure 1) 


y=cssinh &;siny 


x =cs cosh cos 7 = ag cos 9, | 


< 


r 


Figure 1. Notation for elliptic cone. 


To second order the flow is irrotational, so that the velocity is the gradient 
of a potential ®. Introduce a normalized perturbation potential ¢ by 
setting @ = U(s+¢) where U is the free-stream speed. ‘Then the equation 
of motion, including all terms whose etfect may be of second order in the 
slender-body approximation, is 

M?(45 d 2d Pay by Pyy)s (3) 


where 6?=.1/?—1. The boundary conditions are that the flow be tan- 
gential to the surface and that, to second order, the perturbation potential 
vanish at the Mach cone 7 =(x? + y?)!?= ps. 

In R. 'T’. Jones’s slender-wing theory (1946) all terms on the right-hand 
side of (3) are neglected, although for slender bodies with thickness the term 


n= T 
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6?¢,, is implicit in the condition far from the body, as shown by Ward (1949). 
Hence in elliptic-conical coordinates the slender-body equation is simply 


b::+¢,,=0. (+) 
The condition of tangential flow at the surface is found to be 
6.=abs(1+¢,) at €=&, (5) 


where ¢, can be neglected in linearized, and hence also slender-body, theory. 
The general solution of these equations, together with its asymptotic be- 
haviour far from the body, is 


¢=absi+C ~ abslog— +C. (6) 


‘The constant C, which may depend upon s, is evaluated by considering the 
asymptotic behaviour. Ward’s theory provides a general connection 
between the term proportional to log r and that independent of r, which for 
conical bodies is that the perturbation potential should be asymptotically 
proportional to 


a 
‘Thus the slender-body solution is found to be 
a] 

b=abs( §+ | log). (7) 


First derivatives in cartesian and elliptic-conical coordinates are related 
on the surface of the cone by 
bcosy d:—asinn 
* sin*y + 6? cosy)’ 
asinn 6: +b cos 4, 
ry 


(8) 


s(a? sin*7 + cos*y) 
— b?) sin cos 7 — abd: 
s(a? sin*n cos*n) 


+ 


The pressure coefficient is given to second order by 
Cy= —2$,— (47+ $7) + Bb: + + + +47), (9) 
of which only the first two terms are required in linearized theory. Hence 


the slender-body approximation to the pressure coefficient at the surface of 
the elliptic cone is found to be 


4 ab 
C,=ab| 2log —-2+ |, 10 
8 p(a+b) a sin*y + 6? cos*y 
The drag coefficient (referred to base area) is the average with respect to 7 
of the pressure coefficient : 


4 
Cp | C, dy ab| 20g | (11) 


‘These results were given by Fraenkel (1952); equivalent results of greater 
complexity were obtained by Kahane & Solarski (1953). 
A2 
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2.2. The not-so-slender solution 

The slender-body solution is now to be refined by iteration. Rather 
than attack at once all the terms neglected on the right-hand side of (3), 
consider first only the single linear term $74... ‘Thus a second approxi- 
mation is sought to the solution of the Prandtl-Glauert equation, which 
may be written for purposes of iteration as 


= (12) 


rat Dyy 


Evaluation of the right-hand side in terms of the previous first approxi- 
mation yields, in elliptic-conical coordinates, 


3 sin?2y — sinh*2& 
= — b*)s sinh 2€ (13) 
(cosh 2€ — cos 27)* 

lhe usual technique of introducing and (€—7)) as independent 


variables permits a particular integral y% to be found directly by integration: 


2€ — cos 27 


2 cos 27 
ys = 1 — b?)s sinh 2¢( —1 (14) 


It can be verified from the general theory (Lighthill 1954, Van Dyke 1956) 
that this particular integral behaves asymptotically in the manner appropriate 
to vanishing perturbation potential at the Mach cone. Hence the con- 
dition far from the body is sausfied provided that the complementary 
function required to restore tangential flow at the surface vanishes at infinity, 
except possibly for a multiple of the first approximation (7). 

The linearized tangency condition has already been satisfied in the first 


approximation, so the correction consisting of the particular integral ys plus 
the complementary function y must satisfy y.+y-=0 at Possible 
complementary functions that vanish far from the body are 
sinh 2€ 
log 2 (cosh 2£ — cos 27) — 2€ 
cosh 2€ — cos 27 

and the combination of these with the first-order solution (7) that satisfies 
the tangency condition is 


3 Be 3. sinh 2& 
—ab{ 1+ log + 5 leg2 (cosh 2€—cos2y)) |. (15) 


From (8) and (9) the surface pressure coefficient is found to be 


+ ab 
C,= abl 2 log | 


+ B*a*b* 15-2 


+ 


3 4 
2 ab“ B(a+b) 


| 
: 
2 
| | 
9 9 
—Zlogz- + sin2y} |, (16) 
By 2 
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where v?=a?sin*) + 6? cos*7. ‘The corresponding drag coefficient is 


Dol 


a? +b + 
E ab ( ab 2) log Bla 


‘The first terms in these expressions are the previous slender-body results, 
and the second terms the not so-slender-corrections. 


Second-order solution neglecting triple products 

Consider next the non-linear terms in the equation of motion (3) that are 
products of perturbation quantities, the triple products in the last group 
being deferred to the next section. ‘Thus the iteration equation is 


Puy by by:) + + (18) 


In the full second-order theory (Van Dyke 1952), and hence also its slender- 
body counterpart, a particular integral vanishing at the Mach cone that 
accounts for all terms on the right-hand side except that involving (y + 1) is 
given in terms of the first approximation by 


on pe sinh 2€ 


cosh — cos 27 


For the remaining term, a particular integral is given, according to the 
general theory (Van Dyke 1956), in terms of the cross-sectional area S(z), 
by the plane wave 


Yo = — = ats, (20) 


2 


alin 


‘The linear term in the tangency condition (5) has been satisfied in the 
previous solution. Hence a complementary function y is required such 
that (since y..=0) 


ab 


+x: =absd, = | —log at €=§£,. (21) 


B(a+b) sin®y + 6? 


The desired result is a combination of the three constituents of the com- 
plementary function in the previous not-so-slender problem : 


4 Be 


4 ) sinh C 
BC Bla +b) -1(1 ~ cosh 2€ — cos 


+ BX 51 log 2 (cosh 2é — cos 2n) — | (22) 


= 


+ 
Bla +6) 
9 

: 
a 
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From (&) and (9) the surface pressure coefficient is found to be 


C,, = ab(2A+ p) + B2ab[Zabr + 3(a? + b2)A— a —b)? + Lab] 


1 


Q2 
p 


ab 


B(a +b) a@* sin*7 + 
‘lhe two rather complicated last terms in the not-so-slender result (16) have 
here disappeared, so that this more complete result is actually easier to 
compute. The drag coefficient is 


ab(2A + 1) + | + 3(a? + b®)A— —b)? + Sab; 


p24 2 


2.4. The effect of triple products 

The triple products appearing as the last terms in (3) are known to give 
contributions of second order for bodies of revolution, although they have 
been negiected in the few existing second-order solutions for thin wings 
Their effect is found by considering the iteration equation 


Pret Pry t (25) 
This is precisely the iteration equation for the second approximation of 
the Janzen—Rayleigh procedure for plane subsonic flow. Moreover, the 
boundary cenditions are essentially the same as in that problem, since the 
normal velocity at the surface must vanish in both cases, and it suffices here, 


as in the Janzen—Rayleigh problem, to require the velocity disturbances to 


vanish at infinity except possibly for a multiple of the first approximation. 
Use can therefore be made of existing treatments of the Janzen—Ravleigh 
problem, of which one of the most elegant is Kaplan’s (1942) method of 
residues. ‘hat method gives the complementary function as well as the 
particular integral directly in terms of quadratures, which is advantageous 


here because, although the previous complementary functions were readily 
found by trial, that to be obtained next could have been guessed only with 
extraordinary insight. 

Kaplan’s method utilizes the complex variable Z in the plane in which 
the body appears transformed into a circle. In the present case, the elliptic 
cone is the unit circle in the Z plane if 


Z=oexp(§+m), o=—. (26) 


The first-order slender-body solution (7) expressed in terms of Z is (following 
Kaplan’s notation) the real part of 


B(a+ 
4 


{(2Z)=abs| | +log (27) 


} 
1 
3 
t 
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and the complex velocity IV, is given by 
abs 
W(Z)=f(Z)= - (28) 


Then Kaplan’s equation (32), without the i term (which corresponds to 
a uniform flow at infinity), with U=1, R,=1, and a factor M? supplied, 
gives for the second-order complex velocity 


(aby 1+o7/ Z | 
=1M? = 
W(Z,Z)=}M (a+ by s| 2 log 
4 | | a 
22 1 | Z+oa 
2Z 1 


Now if the second-order increment in ¢ is the real part of f,(Z, Z), then 
of, OZ =W, and of, Hence integration gives 


(aby 
=-—1M? 
(a+bp* 
+ 1 Z L+oa 
+log 
E 4Z (+5 
1 +o? Z* 
= y (30 


The constant of integration has been chosen so that this behaves far from 
the body like the first-order solution (27). ‘The last term inside the bracket 
is the particular integral, and the remainder is the complementary function, 
which is very complicated when written in terms of real variables. Some 
computation yields the following increment in surface pressure due to triple 
products which is to be added to (23): 


5 
AC, = —log > 
pv 


a*h? az \ a—v? a a 
> cos '- + — ———— cosh !- — log; 
v 


v?(a® — y a 


2_ 2 2h2(g2 — 
sin 2) : ) 

2 


2v® 


b 


a 
where again v?=a?sin*7 + 6?cos*7. The integrals involved in calculating 
the drag aro treated by contour integration of (Z*— «?) 'log(1+oZ) (1 --oZ) 
around the unit circle, together with integration by parts. ‘Thus the incre- 
ment in drag due to triple products, which is to be ee to (24), is found to be 


3 a2 + 5? 4 
ACp= M?a"b? E log Bath * (32) 


} ‘ 
i 
E 
Fa 
= 
; 
‘ 
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2.5. Reliability and accuracy of the solution 

Setting 4 equal to a yields the second-order slender-body solution for a 
circular cone, which agrees with Br uerick’s result (1949). 

As a check, an independent second-order solution (based upon the 
Prandtl-Glauert equation without the slender body-approximation) has 
been carricd out for a slightly eccentric elliptic cone by perturbing the 
solution (Van Dyke 1952) for a circular cone. ‘Three terms were retained 
in an expansion in powers of the eccentricity, and the resulting surface 
pressure and drag, when expanded into slender-body form, agree completely 
with the expansion of the above solution in powers of the eccentricity. This 
check, which is so elaborate as to be almost conclusive, was applied also to 
the not-so-slender solution. As an additional check, the pressure has been 
studied in the vicini.y of the leading edge as the elliptic cone collapses to a 
flat wing. In this limit the leading edge approaches a thin highly yawed 
parabolic cylinder, and the above solution is found to reproduce the second- 
order Janzen—Rayleigh solution (Imai 1952) for a parabola in the subsonic 
flow corresponding to the component of free-stream velocity normal to the 
edge. 

Despite occasional statements to the contrary, slender-body theory is 
generally somewhat less accurate than linearized theory (of which it is a 
further simplification), particularly at high supersonic speeds, and the 
inferiority is compounded in higher approximations (Van Dyke 1952). 
However, at ‘ordinary’ supersonic Mach numbers, of which the arche- 
type is .V/= \ 2, the first two approximations are known to approach rapidly 
the exact solution for reasonably slender cones (Broderick 1949), and com- 
parable accuracy is to be expected for elliptic cones subtending similar solid 
angles. 

Recently, Rogers & Berry (1956) have measured surface pressures at a 
Mach number of 1-41 over two flat elliptic cones having a= tan 30° =0-577 
and 4=0-05 and 0-10. Figures 2 and 3 show their measurements at zero 
angle of attack, which agree well with the present theory. Pressures are 
plotted against 7 =cos '(v as) in order to expand the narrow region of high 
pressure near the leading edge. Also shown for the thinner wing in figure 2 
is the result of neglecting triple products, which is the second-order theory 
used by Rogers & Berry. Inclusion of triple products is seen to yield 
significant improvement in a small region behind the leading edge. The 
maximum increment in pressure coefficient due to triple products is 0-036, 
whereas that due to double products is — 0-058. For the thicker wing the 
corresponding maximum, values are actually smaller, being 0-027 and 
— 0-032, but the maximum net effect of both kinds of non-linear terms is 
somewhat greater than that for the thinner wing. 

These two wings are equivalent in cross-sectional area to circular cones 
having semi-vertex angles of 9-6 and 13-5 , for which second-order slender- 
body theory predicts pressure coefficients only a few per cent too high 
(Broderick 1949). However, it might be feared that, with Sa=0-577, the 
planform is so wide that the error in the slender-body expansion would be 


if 
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much greater. At least in the limit of vanishing thickness it can be shown 
. that this fear is unwarranted. The surface pressure almost everywhere is 
given exactly in the limit by Squire’s (1947) linearized solution, and second- 
order slender-body theory is 5-7% low. At the leading edge the 
pressure rises to a maximum whose exact value in the limit of zero thickness 
corresponds to loss of the normal component of free-stream velocity, and 
second-order slender-body theory is there only 1-8°, low. (The 
corresponding defects in first-order slender-body theory are 18-8°,, and 
15-2°,.) Thus the good agreement between theory and experiment 


shown in figures 2 and 3 might reasonably have been expected. 


0.5 
0.4 upper surface 
‘al 
x lower surface 
0.1 
\ ane 
2nd-order slender-body theory 
_2nd-order slender-body theory 
without triple products 
et ° x 
2 
Slender- body theory 
0) 10 20 30 40 50 60 70 80 90 
7, degrees 
Figure 2. Theoretical and experimental pressures over 10°%, thick elliptic cone. 
0.5 
© upper surface IN 
~ x lower surface 
0.4 
0.3 
0.2 
0.1 
0 10 20 30 40 50 60 70 80 90 
7, degrees 
Figure 3. Theoretical and experimental pressures over 20%, thick elliptic cone. 
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3. A CRITICAL DISCUSSION OF OTHER APPROXIMATIONS 
3.1. Slender-body and linearized theory 

The above results are now to be used as a model of the full solution, so 
that other approximations will be considered as additional to those already 
made. The second-order slender-body solution will therefore be referred 
to as the ‘exact’ solution. On this basis, the not-so-slender approximation 
is the ‘ linearized ’ approximation. 

In this model, as in the full theory, expansion of the linearized solution 
in powers (and logarithms) of the body thickness yields the slender-body 
solution as the leading terms. Slender-body, linearized, and exact theory 
have frequently been compared for circular cones (e.g., Broderick 1949). 
Consider, therefore, the other extreme of flat cones, as exemplified by the 
wing of figure 3. For convenience, all comparisons are made at W= , 2. 
The surface pressures predicted by slender-body, ‘linearized’, and ‘exact’ 
theory are compared in the left half of figure+. Within this model, linearize’ 
theory is seen to be definitely superior to slender-body theory. (It becomes 
even more accurate if one uses the linearized velocity components in the 
full isentropic pressure-velocity relation.) 


3.2. First-order thin-wing theory 

Thin-wing theory is a simplification of linearized theory in which the 
condition of tangential flow is imposed not at the surface but at (say) the 
plane y=0. Again, in the second and higher approximations the tangency 
condition is transferred to that plane by ‘Taylor series expansion. — It is well 
known that, as a consequence of this planar approximation, the solution is 
not uniformly valid, in particular along the edges of the wing, where singu- 
larities arise that are intensified in higher approximations. 

‘The thin-wing solution for an unyawed elliptic cone lying inside the 
Mach cone in supersonic How was given by Squire (1947), who found the 
surface pressure to be constant. ‘lhe value involves complete elliptic 
integrals, but two terms of their series expansions yield the model required 


+ + 
C. 2ab (log 1) } log 3a -2). (33) 


‘The same result is, of course, obtained by expanding formally the not-so- 
slender solution (16) for small thickness 6 and retaining only linear terms. 
‘This two-term epproximation is within three per cent of Squire’s result 
for semi-apex angles as great as half the Mach angle, whereas the first term 
alone (the thin-wing limit of slender-body theory) gives that accuracy only 
out to one-fifth of the Mach angle. Comparison with the ‘exact’ pressure 
distribution is shown on the right half of figure 4. 

Integration of this pressure over the surface would give a drag coefficient 
equal to the right-hand side of (33). However, R. T. Jones has pointed out 
(1950) that this procedure misses a term associated with the singularity at 


here: 
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the edge. His full expression for the additional leading-edge drag of a flat 


elliptic cone, and its representation in the present model, are = 
ab 
DUE. \ (1— B?a?) 5B ( ) 
0.5 | 
0.4K 


2nd-order 


/ 
J 


0.3 


‘Exact’ (2nd-order | siender-body) 


Linearized’ 


“Ist-order thin-wing' 


or Slender-body 
| 
0 45 90 ISS 180 


n, degrees 
Figure +. Effect of further approximations upon second-order slender-body solution 


for clliptic cone of figure 3. 


Adding this to (33) gives 

abl 2log 1) log — — 5) (35) 

‘This agrees with the result of expanding the ‘ exact’ drag (24 + 32) for small 
b. However, expanding the ‘linearized’ drag (17) gives — 1 instead of — 3 
in the last term. ‘This discrepancy indicates that linearized theory (even 
with tangency imposed at the actual surface) fails to predict the leading-edge 
drag correctly. Non-linear terms are required, which are implicit in Jones’s 
expression (34), although he derived it by an ingenious use of linearized 
theory. 


3.3. Second-order thin-wing theory 
Formal expansion of the ‘exact’ pressure (23 +31) for small thickness 
and retention of squares of 4 yields the ‘ second-order thin-wing’ result 


los 1) + log -2) 
pba 


? 2 
+2- 2log —log ——— +(y+ (36) 
Pa 


Basin 7 B 


Fenain & Germain (unpublished) have recently calculated the full second- 
order thin wing-solution for the elliptic cone. When reduced to the present 
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approximation, their expression for pressure agrees with (36) except for the 
terms ./*a*b*[ (siny) *—log(2 Basiny)] arising from the triple products, 
which they do not consider. 

The contribution of the triple products is of the same order as the other 
second-order terms, and, as pointed out above, the numerical values are of 
the same magnitude. Hence it appears that triple products must be 
included in any complete second-order thin-wing solution for a wing having 
round subsonic edges. ‘ 

Several of the second-order terms in (36) are infinite at the leading edge 
(7, =), so that the pressure distribution has the form shown on the right in 
figure +. ‘The singularity is non-integrable, and must therefore be elimi- 
nated betore the pressure can be integrated to find the drag. Corrections 
of this sort are discussed below. 

Moore (1950) has given the full second-order thin-wing solution for a non- 
lifting cone of diamond section lying inside the Mach cone. For narrow 
planforms and very small thickness he finds large second-order effects over the 
entire wing surface. Lighthill has suggested’ (1954) that this may indicate 
complete breakdown of the planar approximation for narrow wings having 
stagnation edges. However, the round edges of the present example are a 
more severe test of the planar approximation, and the moderate magnitude of 
second-order etfects shown in figure +(except for the inevitable leading-edge 
singularity) suggests rather that Moore’s computations are in error. ‘This 
conclusion is confirmed by Fenain & Germain’s recent recalculation (1955) 
of Moore’s problem, which has uncovered an error in his analysis. 

The ‘ first-order’ pressure coefficient (33) has the form 58 'F,(8a, Bd, 7) 
consistent with the generalized Prandtl-Glauert rule, and the ‘ second- 
order’ increment in (36) has the form ?F,+(y+1)M's 
found by Fenain & Germain(1955) in their treatment of the flat diamond cone. 


3.4. Edge corrections 

Recently the author has proposed rules for rendering thin-wing theory 
uniformly valid at subsonic edges (Van Dyke 1954). The rules for round- 
nosed airfoils are based, to second order, upon consideration of subsonic 
How past a parabola having the same radius as the edge, and consist essen- 
tially in multiplying the thin-wing solution by the ratio of the exact velocity 
on the parabola to its thin-wing expansion. It was suggested that for swept 
edges the rules are to be applied to the normal component of velocity. 

In the present model, the ‘exact’ solution for a parabola in subsonic 
flow consists of the first two terms of the Janzen—Rayleigh approximation 
(Imai 1952). Applying the resulting correction to the ‘first-order thin- 
wing’ solution (33) yields the uniformly valid first approximation 
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Here ‘C,,’ is the formal thin-wing result (33); and the upper and lower 
signs apply to the edges where »=0 and z= respectively. This approxi- 
mation can scarcely be distinguished from the ‘exact’ solution to the scale 
of figure +; the comparison is therefore shown in the following table: 


5 10 20 90 
+380 | 0-317 | 0-220 | 0-145 
)-370 | 0-305 | 0-212 | 0-132 


n, degrees 0 
* Exact ’ Cy, (23-+-31) 0-404 | ( 
Uniform 1st approx. (37) | 0-403 | 0- 


The refinement of including etfects of both edges, not made here, would 
improve the agreement, giving, for example, 0-138 instead of 0-132 at 7» =90 . 

An unexpected complication arises when the corresponding correction 
is undertaken for the ‘second-order thin-wing’ solution (36). Although 
the algebraic singularity is removed, the logarithmic singularity persists in 
the component of velocity tangential to the leading edge, and hence also in 
the pressure. ‘The result is therefore not uniformly valid. This means 
that application of the two-dimensional correction to the normal component 
of velocity at a swept edge, which was thought (Van Dyke 1954) to be 
evidently valid, is incorrect for round edges in the second approximation 
except for incompressible flow. A re-examination of this problem is 
necessary. 


3.5, Linear perturbation of flow past circular cone 

Ferri and his colleagues have made considerable use of the * linearized 
characteristics method’, which consists in linearizing in the departure from 
some known basic flow other than a uniform stream. ‘Thus non-circular 
cones are treated by linear perturbation of the known solution for circular 
cones. Although this approximation implies that the cross-section deviates 
only slightly from a circle, Ferri (1951) has reported good agreement with 
experiment for an elliptic cone of 3: 1 axis ratio. 

‘The present model has been treated by following Ferri’s procedure of 
expanding the velocity components in Fourier series in the polar angle 0, 
linearizing consistently with respect to deviation of the cross-section from 
circular, but then calculating pressure from the full relation, which is (9) 
here. ‘The left half of figure 5 compares the resulting pressures with the 
‘exact’ solution for an elliptic cone of 3:1 axis ratio and area equivalent to 
a 10° circular cone, at. @=, 2. The rather large discrepancy can be under- 
stood by cynsidering the corresponding expansions in the still simpler case 
of slender-body theory. Since tan» =(a@ b)tan 0, the slender-body pressure 
(10) can be shown to have the Fourier expansion in 4 
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Termination of this series at cos 10@ introduces a considerable error, as 
shown on the right half of figure 5. An additional error of roughly the 
same magnitude is seen to arise from linearization in the departure from the 


0.2 


Slender-body to cos 1097 


‘Exact’ (2nd-order siender-body) 


of circular cone Slender-body: linear perturbation 

cos 108 of circular cone to cos !0@ 

305 
0 
O 45 90 iso 180 


7, degrees 
Figure 5. Accuracy of method of linear perturbation of flow past circular cone. 


solution for a circular cone (which in this case affects only the surface boun- 
dary condition, since the equation of motion is already linear). ‘The latter 
error might be reduced by retaining non-linear terms in the boundary con- 
dition, as has been suggested (Ferri, Ness & Kaplita 1953). However, the 
remaining error inherent in curtailing the Fourier series is so great that if 
only terms up to cos 104 are to be retained, it would seem that the method 
must be restricted to more nearly circular bodies. 


‘This work was carried out at Cambridge University while the author 
held a Fulbright research grant and a John Simon Guggenheim Memorial 
Fellowship. He is indebted to E. W. E. Rogers for helpful discussion and 
for making his and Berry’s experiments available. 
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SUMMARY 
‘This paper proposes a theory of collisions between small drops : 
in a turbulent fluid which takes into account collisions between é 
equal drops. ‘The drops considered are much smaller than the i 
small eddies of the turbulence and so the collision rates depend i 
only on the dimensions of the drops, the rate of energy dissipation « f 
and the kinematic viscosity v. Reasons are given for believing 
that the collision eficiency for nearly equal drops is unity, and the 
collision rate due to the spatial variations of turbulent velocity is 
shown to be \V= 1-30(r, v)'*, valid for +, r, between 
: one and two. A numerical integration has been performed using 
this expression to show how an initially uniform distribution will i 
change because of collisions. An approximate calculation is then : 
made to take account also of collisions which occur between drops 
of different inertia because of the action of gravity and the turbulent k 
accelerations. 


f The results are applied to the case of small drops in atmospheric 
clouds to test the importance of turbulence in initiating rainfall. 
Estimates of « are made for typical conditions and these are used to 
calculate the initial rates of colhsion, the change in mean properties 
and the rate of production of large drops. It 1s concluded that i 
the effects of turbulence in clouds of the layer type should be small, 
but that moderate amounts of turbulence in cumulus clouds could 
be effective in broadening the drop size distribution in nearly 
uniform clouds where only the spatial variations of velocity are 
important. In heterogeneous clouds the collision rates are 
increased, and the effects due to the inertia of the drop soon 
become predominant. ‘The effect of turbulence in causing 
collisions between unequal drops becomes comparable with 
that of gravity when is about 2000 cm? sec *. 


1. INTRODUCTION 
It has been agreed for many years that while the initial process in the 
formation of clouds i the atmosphere must be one of condensation from 


the vapour phase, this process is not sufficiently rapid for the small water 
droplets to grow to raindrop sizes in the times usually available. Bergeron 
(1933) suggested that ice crystals might play a crucial part in the mechanism 
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by which raindrops are formed, but observational evidence since that time 
has proved that heavy rain can fall from clouds whose temperatures are 
nowhere below the freezing point. 

Various workers have therefore developed the idea of the growth of 
raindrops by the coalescence of liquid cloud drops. ‘These theories suppose 
that drops larger than the mean drop size fall through the cloud under the 
action of gravity, and sweep up some of the smaller drops in their path. 
The rate of growth of the falling drops can be calculated and is found to be 
a rapidly increasing function of the size of the drop. 

Squires (1952) has investigated theoretically the size distributién of 
cloud drops condensing on hygroscopic nuclei, and has concluded that 
only under rather special circumstances can a few droplets considerably 
larger than the mean be formed. If, for example, there are a few giant 
nuclei, these will quickly grow, by condensation, to such a size that the 
above coalescence mechanism will become important. Otherwise, the 
calculated drop size distribution is more uniform than that observed, as is 
also clear from the work of Howell (1949), and the mean size is smaller than 
that for which the above coalescence mechanism becomes more important 
than condensation. 

A natural suggestion therefore is that turbulence in a cloud might lead 
to collisions among the drops, thus giving a greater spread of cloud drop 
sizes and providing the larger droplets on which raindrops could form. 
(Bowen (1950) has indeed considered the history of drops of double the 
mass of the rest of the drops, these drops being supposed formed by the 
random collision of two drops, but the mechanism of such collisions was not 
examined.) East & Marshall (1954) have discussed the previous (mostly 
qualitative) suggestions which have been made about the role that turbulence 
could play in prec‘pitation, and they propose a new theory in which the 
effect of random motions is regarded as being equivalent to the action of an 
increased gravitational field. It is concluded that turbulence could be 
important in a heterogeneous cloud if the random air acceleration is com- 
parable with the acceleration due to gravity. It should be noted that the 
process pictured by these authors is still ineffective if the drops are small, 
and that they predict zero collision rates in homogeneous clouds. 

It is the purpose of the present investigation to discuss a mechanism 
due to turbulence which gives collisions between equal drops and which 
employs a more realistic model of the nature of the turbulent motion than 
that used by East & Marshall. The collision frequencies are found in 
terms of the rate of energy dissipation per unit mass in the cloud and this 
quantity is estimated from the large scale properties of the air motion inside 
the clouds. 


2. 'THE NATURE OF THE TURBULENT MOTION 


East & Marshall regarded the turbulence as equivalent in its effects to 
a random motion in time of the whole air parcel containing all the drops, 
and neglected the spatial variations which are surely an essential feature of 
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turbulent motion. All collisions on their model are therefore due to the 
different motion, relative to the air, of drops of different sizes; and hence 
they found that equal drops did not collide. The effect of the spatial 
variations, however, is to give neighbouring drops different velocities and 
thus cause collisions, whatever the ratio of the sizes of the drops. 

The cloud droplets that we are considering are usually smaller by at 
least an order of magnitude than the length scale of the small eddies of the 
turbulence, and so the relative motion of two neighbouring drops will be 
governed by the small scale motion. It has been pointed out by Batchelor 
(1950) and others that because the Reynolds number for turbulent motion 
in the atinosphere is usually large, the similarity theory of turbulent motion 
will hold for the smal! scale motion. This theory implies that for scales 
of motion sufficiently small compared to the energy-containing eddies, the 
motion is isotropic and the mean values of quantities related to the turbulence 
will depend only on the kinematic viscosity v and the rate of energy dissipation 
per unit mass «, provided that the quantities concerned depend strongly on 
the small scale properties of the turbulence. For example, the relative 
diffusion of a cloud of smoke can be treated in this way, whereas diffusion 
from a fixed source can not, since in the latter case the large eddies are also 
important. It follows that the effect of turbulence in causing collisions 
between neighbouring droplets will also depend on the rate of turbulent 
energy dissipation per unit mass e, and it will be necessary to make an 
estimate of this quantity in clouds under various conditions. We shall do 
this first before describing the mechanism of collisions. 

runt (1939) has given as an average value in the lowest few kilometres 
of the atmosphere «=5cm?sec *. Few direct measurements which allow 
« to be estimated in clouds have been made, but R. J. ‘Taylor (1952) has 
deduced values of « of the order of 1000cm?sec * close to the ground in 
moderate winds. Another approach to the problem is available however. 

The results of laboratory experiments (e.g. Batchelor & ‘Townsend 
1948) have indicated that the rate of turbulent energy dissipation is usually 
of order u?/l, where u is a root-mean-square turbulent velocity and / is a 
length scale associated with the energy-containing eddies. Measurements 
of accelerations experienced by aircraft in bumpy conditions indicate that 
fluctuating velocities of a few metres per second are common, with corre- 
sponding eddy sizes of the order of tens or hundreds of metres. If we take 
as typical figures, = 2 metres sec, /= 50 metres, we obtain «= 1600 cm? sec 
We might suppose that something less than the mean figure of « = 5 cm?sec 
would be applicable to stratiform clouds where there is a small mean 
velocity, and take a larger value say « = 1000 cm?sec * as an estimate of the 
conditions in turbulent cumulus clouds. However, it may be possible to 
make better estimates than these for the different types of clouds. 


3. DiIscUSSION OF THE COLLISION PROCESS 


‘There are two ways in which turbulence causes collisions between 
neighbouring droplets. Firstly, there are the spatial variations of the 
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turbulent motion referred to previously. Collisions due to this process 
can conveniently be called ‘collisions due to the motion of the droplets 
with the air’. Secondly, each droplet is moving relative to the air sur- 
rounding it, owing to the fact that the inertia of a droplet is different from that 
of an equal volume of air. It follows that neighbouring droplets of unequal 
size will have different velocities (since the inertia of a droplet depends 
on its size), and this also will lead to collisions. ‘This process is called 
‘collisions due to the motion relative to the air’. ‘The latter process will 
not give collisions between droplets of equal sizes, such collisions being 
due to the first process. 

Before we can calculate the rate at which collisions between droplets 
occur, it is necessary to consider the effect that the presence of a droplet 
has on the motion of neighbouring droplets, that is, we must consider the 
distortion of the flow due to the presence of a drop. A measure of this 
distortion is the collision efficiency, which can be defined as that proportion 
of drops which would have collided in the absence of distortion, actually 
to do so. It 1s clear that the collision efficiency must be dependent on the 
nature of the flow. Collision efficiencies were calculated by Langmuir 
(1948) for the case of small drops suspended in the steady laminar flow 
around a fixed large sphere. His numerical results are only likely to be 
accurate in cases where the basic assumptions are satistied, although they 
have been extensively applied beyond their range of validity by others, and 
by East & Marshall in particular. It is difficult to justify the application 
of Langmuir’s results to the present problem, since we are concerned with 
a case in which the colliding drops are of nearly equal sizes. Recent 
measurements, by ‘Telford, ‘Vhorndike & Bowen (1955), of collection 
efficiencies for nearly equal drops somewhat larger than cloud drop sizes, 
have indeed given values considerably higher than would be predicted by 
an application of Langmuir’s theory. 

Now the Reynolds number of the relative motion of two nearly equal 
approaching cloud droplets will usually be much less than one, and some 
relevant experimental evidence for this range of Reynolds numbers is 
provided by the experiments of Manley & Mason (1952, 1955). ‘They 
observed the collisions between glass spheres and between air bubbles sus- 
pended in a uniformly sheared viscous liquid, and showed that in these 
circumstances, the collision efficiency is indeed unity, that is, the distortion of 
the flow does not influence the collision rate. 

It thus seems that, in the absence of further evidence, it is not unreason- 
able to take the collision efficiency of nearly equal droplets as unity; this is 
equivalent to neglecting altogether the distortion of the flow by a drop. 

We are interested primarily in collisions between nearly equal drops, 
and as a beginning we shall confine our attention to collisions due to motion 
withthe air. Collisions due to motion relative to the air are not unimportant, 
but this process does not give collisions between equal drops. Later, we 
shall mzke an estimate of the relative importance of the two processes. 

It is useful to keep in mind the qualitative picture of the local shearing 
motions which lead to collisions between drops carried along with the fluid. 
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; For two close points in a turbulent fiuid, the relative motion is that of 
: uniform strain. ‘Taking an origin at the centre of one drop and the 
. coordinate axes in the directions of the principal rates of strain, the stream- 
lines of the relative motion in one of the coordinate planes are as shown 
in figure 1. The other drops are moving with the air along these stream- 
lines and since we neglect the distortion of the flow field due to the presence 
of the drops, for the reasons mentioned above, the collision rate of the 
‘fixed drop’ with other drops is just the flux of fluid inwards across the 
; surface of a sphere, concentric with the fixed drop and of radius equal to the 
a sum of the radii of the two approaching drops, multiplied by the number 
density of the other drops. ‘The equivalent calculation for a uniform 
laminar flow assumes, of course, that all drops in a cylinder parallel to the 
flow and containing the effective cross-section of the fixed drop should meet 
ra that drop. 


Figure 1. Streamlines of the relative motion in one of the principal planes. 


+. COLLISIONS BETWEEN DROPS MOVING WITH THE AIR 


We now proceed to the calculation of this collision rate. Let the mean 
concentrations of two sizes of drop in a given population be 2, and n, per 
unit volume, and their radii 7, and 7, respectively. ‘The collision radius 
for a pair of drops, one of each type, will be just the sum of the two radii, 
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R=r,+1r, say. Now the mean flux of fluid into a sphere of radius R 
surrounding one drop is 


where w, is the radial component of the relative velocity, and the bar denotes 
a mean over many realizations of the motion. If we suppose the turbulence 
to be isotropic, this mean will be equal to the mean at a fixed point in space, 
since there is no correlation between the position of a drop and the properties 
of the turbulence. If the drops are randomly distributed and moving with 
the air, the collision rate is thus 


— | 


w,dS, 
wp <0 
and, by the similarity hypothesis, the integral will be a function of € and v, 
and of the length R. 
To evaluate the integral, we note that the equation of continuity of the 


fluid gives 


w,dS + | w,dS=(. 
Hence, 
entire sphere entire sphere 


and since the small eddies are isotropic, the last expression becomes 
2nR*| w,, where w, is the radial relative velocity along the radius parallel 
to the x-axis. Since RF is usually small compared with the length scale of the 
small eddies (as will be shown in §6), |w,]=R|éu/dx|, where u is the 
x-component of the velocity. 

Further, the mean square of the velocity gradient is related to « and v 
through the expression (du dx)?=e 15v (Taylor 1935). We now assume 
that du 6x is normally distributed and obtain |@u/Ax|=(2e/15v)"".  (du/dx is 
only approximately normally distributed but the error so introduced in 
the present context is not large; see Townsend (1947).) The collision 
rate is thus 


\ 1/2 1/2 


15v 


= Say. 2 (1) 


We can now use (1) to find the rate of growth of a population of equal 
drops. Let us suppose that initially we have a uniform size distribution of 
drops of type 1. Larger drops will form by multiple collisions, and we 
denote by n, the number of drops with s times the mass of a type 1 drop. 
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We also suppose that when two drops collide, they coalesce. Then the 
following equations hold (Smoluchowski 1917), 


an, » 

dt 

dns 

—— = — No — — 
dt 

at 

in 


Equations (2) can be solved in closed form if x;; = 2z, that is, if it is assumed 
that the collision rate is independent of the particle size, and the solution is 
n, No 


(3 
(1+ xnyt)* ! 


where is the concentration (total number of drops per unit volume) after 
time ¢, m) being the initial value. . This solution was employed by Schumann 
(1940) in his study of the size distribution of fog particles, although he did 
not there discuss the effects of turbulence quantitatively. However, it is 
not satisfactory to neglect the dependence of collision rate upon drop size 
without further investigation and we therefore integrated equations (2) 
numerically. It was found that (3) can indeed lead to a serious under- 
estimate of the number of large drops formed in a given time. 

The equations (2) may be rewritten in a convenient dimensionless form 
if nym, and K=1-30tri(e v)!2n, are used as variables. ‘Thus all the 
quantities which are given in a particular application are absorbed into the 
quantity K. The numerical integration was carried out over short intervals 
of K=6~ 10 * by an iterative process, assuming a value for mean concen- 
trations during the interval, working out the number of drops of each size 
produced and removed, and checking the assumed values. 

The integration was continued until the total number of drops was 
reduced to two-thirds of the initial value, that is, until the mean drop mass 
increased by 50°,. At this stage, K 1-044 ‘The effect of multiple 
collisions in broadening the drop size distribution is shown in figure 2, 
where the results of our numerical integration are plotted. ‘Ihe numbers 
of drops of each size up to eight times the original mass are plotted on a 
logarithmic scale for three values of AK. For comparison, the numbers 
calculated assuming a collision rate independent of size (i.e. using (3)), are 
also shown. (The lines in the figure are intended for guidance, and have 
no meaning in themselves.) It is clear that this approximation under- 
estimates the number of eightfold drops present at a given time by a factor 
of twenty or more, and that the discrepancy increases with size. As a 
method for estimating the rate of change of the mean properties, the 
approximation is reasonable provided K is not too large. 
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Let us now use the result of the numerical integration to estimate the 
time necessary for the mean mass of cloud drops to increase by 50°,. 
have 1-30tr3(e v)!2n, = 1-044 x 10 in c.g.s. units, or wt(e v)!? = 3-36 x 10°, 
where w is the liquid water content in gm cubic metre and ¢ is the time in 
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10°? 


Figure 2. Distribution of drop sizes produced from an initially uniform population. 
Circles show the result of the numerical integration, and crosses show the 
result of a calculation assuming the collision rate to be independent of size. 
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seconds. So with typical figures of w=1-5gm cubic metre and v=0-17 
cm? sec |, we obtain as the time for the mean mass to increase by 50°, owing 
to the assumed motion with the air, 
4-() x 104seconds= 11 hours, if «=5cm?sec 
t=2-8 x 10%seconds= hour, if «=1000cm*sec *. 


‘Thus, the motion of the drops with the air will affect the mean properties 
of a cloud rather slowly, except under conditions of vigorous turbulence. 

In applying these ideas to the initiation of rainfall, however, we are 
interested primarily not in how the average drop size changes, but in the rate 
at which a few large drops form. Let us, for example, see how frequently 
the double drops required by the Bowen theory will form. He assumed 
uniform drops of 10; radius, and a liquid water content of | gm cubic 
metre, corresponding to a concentration of 240 drops cm’. ‘The rate of 
production of double drops initially is then 3-0 10 “e cm “sec 
This shows that with « =5cm?*sec °, 0-04°, of the total number are double 
drops after one minute, and with «=1000cm?sec%, 0-6°., are double 
drops after one minute. ‘Thus moderate amounts of turbulence should 
produce a significant concentration of double drops. 

Again, if we apply the numerical results for multiple collisions to a case 
in which =600 drops cm* initially, «= 1-5 gm cubic metre and «= 100 
cm*sec °, we find that in 16 minutes there will be 100 drops per litre with 
mass equal to or greater than four times the original mass, and after 90 
minutes, 100 drops per litre with mass equal to or greater than eight times 
the original. For sizes greater than this, the rates of growth increase 
rapidly, but so too does the relative importance of other processes, as will 
be seen in the following section. 


5. MOTION RELATIVE TO THE AIR 

In this section we shall consider approximately the motion of the drops 
both with the air and relative to the air, and also the relative motion of the 
drops due to gravity. In order to do this, we shall use a model of the 
collision process which will enable us to take all three processes into account 
at the same time. ‘This will provide us with an indication of how the effects 
of the three processes should be added in a case in which they are calculated 
separately in some other way. It is to be noted, however, that the assumption 
that the collection efficiency is unity (see §3) is only likely to be valid for 
nearly equal drops, and so the following calculation will not applv to drops 
of very different sizes*. 

The argument will proceed as follows. The probability distribution of 
the relative velocity of two drops is in general a function of their separation. 
When the separation of the centres of the drops is small compared with the 


*'The experiments of Manley & Mason indicated a change in the relative motion 
of two colliding drops when the ratio of the larger radius to the smaller was greater 
than two. When the ratio was less than two, they observed the collision efficiency 


to be unity. 
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length scale of the small eddies (it is shown in §6 that in clouds the drops 
are usually small enough for this to be possible), this probability distribution 
can be deduced from the equation of motion of a drop and a knowledge of 
the statistical properties of the small eddies of the turbulence, provided 
we neglect the distortion of the flow by a drop for the reasons given in § 3. 
We shall then have, in particular, the probability distribution of the relative 
velocity of two drops just before they collide. We denote by P(w)dw 
the probability that the relative velocity w of two drops just before they 
collide lies in the range w, w+ dw. 

Let us calculate the rate of collision between drops of radius 7, and 
ry, Whose number densities are 7, and x, respectively, in terms of P(w). 
A collision will take place when the separation cf the centres of the drops 
is R=r,+r,. We shall consider the collisions that occur in a small interval 
of time 5/, where df is small compared to the time scale of the small eddies, 
and we may suppose that the drops move in straight lines with constant 
velocity during the interval 6/.__ For the purpose of the calculation, we may 
suppose one of the colliding drops, of type | say, to be at rest, the velocity 
of the other drop just before the collision being w. Then a collision will 
take place in the interval df if, at the beginning of the interval, the centre of 
the second drop was somewhere in a volume 7R?wét, where w=|w]|. We 
now integrate over all values of w and multiply by the number density of 
drops of type 2 and thus obtain the number of collisions between the 
‘stationary’ drop and the other drops in the time of. On multiplying by 
the number density of drops of tvpe 1 and dividing by 45t, we obtain the 
collision rate as 


N=7Rn,n,| | | P(w) dw. (+)* 


It is permissible to multiply by the number density in the above manner 
only when the mean velocity of two colliding droplets is statistically inde- 
pendent of their relative velocity ; for if this were not so, the local concentra- 
tion of droplets would be related to the relative velocity. However, the 
mean velocity is controlled by the large energy-containing eddies and the 
relative velocity by the small eddies, and it is known that the large and small 
eddies are statistically independent when the Reynolds number of the 
turbulence is large (see ‘Taylor 1938). 

We see from equation (4) that the collision rate can be found as soon as 
P(w) is known. P(w) contains the effects of the spatial variations of 
velocity in the fluid (giving collisions due to the motion with the air), the 
turbulent accelerations (giving collisions due to the motion relative to the 
air), and gravity. It could be expressed formally as a complicated function 
of these three effects, but in that form it would not be possible to evaluate 
the integral in (4) in finite or simple terms. Now the variance of w can be 
evaluated without difficulty from the equations of motion of the drop and 

* Further details of the argument leading to (4) may be found in Chapman & 
Cowling, Mathematical Theory of Non-Uniform Gases, Oxford University Press, 1952, 
pp. 60, 89, 90. ‘ 
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the statistical properties of the turbulence; and the method we shall use in 
order to avoid the above difficulty and to give integrals which we can 
evaluate, is to replace P(w) by a simpler function which has the same 
variance. We shall justify this method by showing that, in the special 
cases when either motion with the air or gravity alone is effective, it does 
not lead to significant error. : 

Before discussing the explicit form that we shall take for P(w), it is 
convenient to calculate the variance of w. ‘To do this, we consider two 
droplets which have velocities ¢, and ¢, just before they collide. We denote 
by u, and u, respectively the ‘ undisturbed’ velocities of the air surrounding 
these drops, so that the relative velocities between the drops and the air 
surrounding them are q,=¢,—U,, Q.=¢,—U,. Further, w=c,—¢. 
It follows that the variance of w is given by 


var (w) = var (c, — ¢,) = var (q. — G,) + var (u, —U,), (5) 


provided we suppose u,—u, to be statistically independent of q, and qg. 
lhe second term ts easily calculated since 


var (Uy — = R? (BY 


in view of the isotropy of the small eddies, 
To calculate the first term, we use the equation of motion of a drop, 
which is 


de =. p du p 
u)+ Po dt + — (6) 


where d dt denotes differentiation following a drop, g is the acceleration 
due to gravity, 7 1s the relaxation time of the droplet (+ = 2r29, 9 for small 
spherical drops obeying Stokes’ law), 7 is the radius of a drop, py is the 
density of the drops, p that of the air, and . the viscosity of the air. Hence, 
dq, 1 p \/du 
-1)(F +8). (7) 
dt 1 dt 
Now if 7 is small compared with the time scale of the smallest eddies 
(and this will be shown in §6 to be usualiy the case), we may neglect the 


first term on the left-hand side of (7) and replace (du dt)? by (Du Dt), 
and obtain 


We have already taken account of the relative motion of drops due to the 
difference between u, and u,, so that in considering q, and q, we can suppose 
that the air velocities near two close particles are the same. (This is 
equivelent in effect to the assumption inherent in the work of East & Marshall 
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when they considered a random motion of the air parcel as a whole.) Thus, 
the correlation coefficient between q, and q, is unity and 


/Du\2 
var -£) EG) 


Equation (5) then gives 


var(w)=3(1- 2) (71-7 +(1-4) 7 + . (8) 


Po \ 


We now return to the discussion of P(w). We shall replace P(w) by 
the function 


P*(w) = (2) exp(— Bw’), (9) 


where f is chosen so that the variance given by this distribution is the 
same as that for P(w), that is, 38 ! is equal to the expression (8). ‘This seems 
to be the function that corresponds most closely to physical reality and vet 
renders the integration of (4) simple. If the relative velocity were due 
only to motion relative to the air, it seems likely that (9) would be accurate. 
Further, although the use of (9) replaces the constant relative velocity that 
would arise from the differential fall of droplets of different sizes under 
gravity alone by a random relative velocity with the same mean square, it 
will be seen below that the error so introduced is not large. 

When we consider the relative velocity due to the spatial variations of 
velocity in the fluid, however, we find that the choice of the distribution (9) 
is not free from inconsistency. ‘l’o see what is happening, let us suppose 
that the droplets are moving entirely with the air, that is, that the relative 
velocity due to gravity and motion relative to the air is zero. It seems that 
in this case (9) is inconsistent with the isotropy of the small eddies. For 
instance, if the separation of the two droplets is parallel to the x-axis, it 
follows from the isotropy that (w,)? = (w,)? = 1(w.)?; but according to 
equation (9), (wz)? =(w,)? = (w,)?. However, despite this anomaly, the 
result obtained by the use of (9) will be shown to be not too different from 
that obtained in $4, where the model used is consistent with the isotropy 
of the small eddies but still assumes the normal distribution of the separate 
components of the relative velocity. 

We now use (9) to calculate the collision rate. On substituting in (3) 
and integrating, we obtain 


N=7R*n,n, | | in(=) exp (— Bw!) dw= 
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‘The first two terms are zero when 7, and 7, are zero and the drops move 
with the air, when 7,=7, and the drops are identical, or when p =p, and the 
drops have the same density as the ambient fluid. When one of these is 
the case, (10) gives 


€\l2 
) = 167 ‘ 


It will be noticed that the constant here differs from that of equation (1), 
for the reasons mentioned above. Further, in the case when there is no 
turbulence, (10) gives 


0 


whereas the exact calculation for this case, wuch is easily done, gives 


N= ( 1 |r, —72]. 
Po 


Again we note that the difference in the constant is quite small. 


6. COMPARISON OF THE VARIOUS COLLISION MECHANISMS 


We now use (10) to compare the relative importance of the two turbulent 
collision mechanisms in clouds. It must first be shown that the radius 7 
and relaxation time += 2r?p, 9 of the droplets are usually small compared 
with the length scale (1? «)!4 and time scale (v €)!* of the small eddies, 
respectively, since the arguments leading up to (10) depend upon these 
conditions being satisfied. Now for drops of radius 74, 7=6:3 x10 #4 
seconds, and for drops of radius 20 7=5-2 x 10 3 seconds. ‘The corre- 
sponding values of z(e v)'? and r(e,+°j'* are shown in the table. It can be 
seen that assumed conditions will be satisfied except perhaps for the largest 
drops in very strong turbulence. 


| 

«=S.cm*sec™ | 23x10" sec 3°8 x 10-4; 1-1 x 10> 
| 

1000 cm? sec 4°93: 10-7 | 1000cm?sec *| 1-4 102 4-0 x 10-* 


‘lhe expression (10) may be evaluated in terms of known quantities by 
means of a result obtained by Batchelor (1951), who showed that 


a) = 
Dt 


approximately when the Reynolds number of the turbulence is large 
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Hence, the ratio of the first term to the third in (10) is 


2 12 
11-7(1 (7 — 79)? = 8-2 x Ar, 
for water drops in air, where 7, and r, are the drop radii in centimetres, the air 
viscosity has been taken as 10 4c.g.s. units, and p as 1-00 x 10 
gm,cm* (values at 0°C and 800mb). ‘The ratio of the two terms is unity, 
that is, the effects are comparable when r,—r,=2-3 1 for «=5cm?sec 3, 
and when 7, —7r,=0-6 for = 1000 cm? sec 

Thus, except for the very small cloud droplets, the collision rates will 
increase to several times the values in §4 when one drop is a few times the 
mass of the other; and eventually the collisions due to turbulence in hetero- 
geneous clouds will be dominated by the motion relative to the air. (It 
should be emphasized again that the theory developed here is not likely to 
remain valid when the droplets become very different in size.) 

It is also clear from (10) that the motion due to gravity and the motion 
of the drops relative to the air due to turbulence become comparable in 
their effects when the root-mean-square acceleration is of the same order 
as the acceleration due to gravity g; using (10) the condition for equality 
of the effects is g=2v "434, which corresponds to «=2100cm?sec *. 
Thus what we have called vigorous turbulence, with «= 1000cm?sec °*, 
would be less effective than gravity in causing collisions. 

Reviewing the results then, we have shown that if the low collection 
efficiencies calculated for uniform streaming and two drops of very different 
sizes are rejected as irrelevant for the case of nearly equal drops, considerable 
rates of coalescence are predicted in turbulent air, even in clouds with very 
uniform size distributions. In this latter case, the first collisions are due 
to the drops moving with the air, a process which does not seem to have 
been taken into account in previous theories. Particularly for the small 
drop sizes, where condensation is also important, the collisions of nearly 
equal drops concurrently with the condensation could lead to a broadening 
of the distribution of cloud drop sizes. As the size increases, and the range 
of sizes becomes larger, collisions due to the different relative motions 
between drops and the air (induced by the turbulence and a steady fall 
under gravity) will become predominant. Our theory does not allow us to 
follow the growth of a drop right up to raindrop sizes. 

The magnitudes estimated for the initial collision rates are such that the 
properties of stratiform clouds would be affected rather slowly, whereas 
in a vigorous cumulus updraught the rate of production of larger drops 
might be sufficient to initiate the formation of rain. ‘Thus we are not faced 
with the embarrassment of predicting that all clouds should fall out as rain. 


The authors wish to express their thanks to Dr G. K. Batchelor for his 
suggestion of this problem, and for many helpful discussions. ‘The work 
was made possible for one of us (P. G. S.) by the award of a Department of 
Scientific and Industrial Research maintenance grant, and for the other 
(J. S. T.) by the tenure of an Exhibition of 1851 Overseas Scholarship. 
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SUMMARY 

Sir Charles Darwin has advocated a study of the ‘drift’ of 
material surfaces in the classically investigated irrotational flows 
past bodies. This suggestion is followed up and given further 
support in the present paper. In particular, it is shown how 
secondary flows can be evaluated by use of the ‘ drift function’ ¢ 
for the primary flow. ‘This is a function ($1) such that material 
surfaces initially at right angles to the stream drift into shapes 
expressible by equations ¢ = constant. 

The analysis leads to a simple expression (§ +) for the secondary 
velocity field in the flow past an infinite cylinder of any cross- 
section, with the upstream velocity normal to its axis and increasing 
linearly with distance along the axis—a problem in which only the 
secondary vorticity field was previously known. ‘The drift past a 
sphere is computed and illustrated (§ 6), and the secondary vorti- 
city field in shear flow past a sphere is tabulated (§7). ‘There is: 
also a detailed study (§ 3) of the asymptotic form of the secondary 
velocity field in flow past any body, based on a result of Darwin 
concerning ‘ hydrodynamic mass’. 


1. INTRODUCTION 

Sir Charles Darwin (1953) has shown what interesting and important 
conclusions can be drawn by studying the deformation, or ‘drift’ as he 
calls it, of material surfaces—-and, generally, by studying the motions of 
individual particles—in the classical problems of irrotational flow of fluid 
about bodies. Some more advantages of this approach are indicated in the 
present paper. In particular, the Cauchy—Helmholtz—Kelvin result that 
vortex lines move with the fluid, while the magnitude of the vorticity 
changes in proportion to the local stretching of the vortex lines, is combined 
with it to deduce how the vorticity in a weakly-sheared flow is altered in the 
presence of an obstacle, as a result of deformation of vortex lines by the 
irrotational component of the flow about the obstacle. 

Detailed results are obtained for the sphere, which Darwin (1953) 
treated only briefly. ‘The deduced vorticity field in shear flow past a 
sphere (§ 7) will be used in a later paper, together with the recent derivation 
(Lighthill 1956) of the image system of a vortex element in a sphere, to 
compute certain features of the secondary flow, and in particular the up- 
stream displacement of the stagnation streamline——a quantity related to the 
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‘displacement of the effective centre’ of pitot tubes in shear flow, des- 
cribed by Young & Maas (1936). 

The approach leads also (§ +) to a convenient expression for the secondary 
velocity field in the flow past an infinite cylinder of any cross-section with 
the upstream velocity normal to its axis and increasing linearly with distance 
along the axis; in this problem, only the secondary vorticity was previously 
known (Hawthorne 1954). 

Hydrodynamics has achieved impressive results by the simplification of 
concentrating on the steady (or nearly steady) field of flow velocities relative 
to a moving body, as specified in the Eulerian manner. ‘This does not 
mean, however, that additional information, regarding the history of 
individual particles of fluid, is of no value. For example, Darwin (1953) 
shows that, when a circular cylinder moves through otherwise undisturbed 
Huid, a typical fluid particle begins to move forwards before the body reaches 
it, then starts to move outwards (away from the path of the body) and 
forwards, then outwards and backwards, and is moving backwards (parallel 
to the path of the body but in the opposite sense) when the body is abreast of 
it; it then moves inwards and backwards, inwards and forwards, and finally 
comes to rest while moving forwards, at a point slightly ahead of its original 
position. ‘This result throws light on the familiar gyrations of the air, and 
of small particles suspended in it, when any large body moves past. 

Darwin considers also “an infinite thin plane of fluid, at right angles to 
the motion, marked so as to be made recognizable, perhaps by means of 
some dye-stutf’’, and asks “‘ after the passage of the body, what is the form 
assumed by this infinite plane ?”’. The answer is that the part near where the 
body has been has drifted forwards, and that the fluid between the initial and 
tinal positions of the material surface has a mass equal to the *hydro- 
dynamic’ or ‘ virtual’ mass associated with the body’s motion. 

It could, perhaps, have been expected that the hydrodynamic mass would 
equal the net ‘mass movement’ (integral of momentum with respect to 
time) induced in the fluid per unit distance travelled by the body, so that the 
total mass movement per unit distance zquals the mass of the body plus 
the hydrodynamic mass. However, the study of hydrodynamic mass from 
an Eulerian point of view had left this interpretation hidden until Darwin 
discovered it. 

It should be emphasized that no abandonment of the Eulerian steady 
velocity field is involved in the studies here described, nor is there any 
question of using the Lagrangian equations of motion. Rather it is neces- 
sary to determine, not only the streamlines, but also on each streamline 
the time at which a fluid particle reaches any given point, measured from 
some fixed time for the particle—say, from when it passes across a particular 
plane at right angles to the undisturbed stream (as in Darwin 1953), or (as 
in the present paper) from when it would have passed across a particular 
plane had the stream remained undisturbed. ‘Thus, if the latter plane is 
y=0, then we require the solution of 
with >Oasx---—x (1) 
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(where wv, is the x-component of velocity, etc., and the undisturbed stream 
has v,=U,v,=v,=0). The advantage of this definition over Darwin’s 
is that, far upstream, material planes at right angles to the stream are planes 
t=constant; hence the surfaces into which these are deformed, in which we 
are particularly interested, are all given by ¢= constant. 

It is easy, in terms of this ‘drift function’ ¢, to determine how vortex 
lines, initially at right angles to the stream, are deformed and stretched by 
the irrotational flow about the body (here, the vortex lines are supposed 
weak enough for the additional stretching by the secondary flow due to 
their presence to be neglected). A vortex element joining a pair of points 
on two neighbouring streamlines will always join points on these two stream- 
lines, and they will always be points with the same value of ¢. ‘The vector 
separation of these two peints, divided by the scalar separation of the 
original pair, gives the direction and magnitude of the vorticity divided by 
its upstream magnitude. The details come out simply once the drift 
function ¢ is known. 

The work is carried out in detail for a sphere, as being the three- 
dimensional shape for which (owing to the simplicity of the image system) 
the secondary flow has the best chance of being calculated. ‘This secondary 
How will be studied in detail in a later paper. Here, only its behaviour at 
large distances from the sphere is worked out, since this deduction does not 
depend on the details of the image system, but conversely leans heavily on 
Sir Charles Darwin’s result concerning hydrodynamic mass. In fact, it is 
shown (§3) that, for any finite body so placed that the irrotational flow 
exerts no yawing moment on it, the part of the secondary velocity field which 
falls off like the inverse square of distance from the body has a simple form 
depending only on the hydrodynamic mass associated with the body’s 
motion and on the mass of fluid displaced by the body. 

Finally, it is pointed out that this secondary velocity field induces what 
might be called a ‘tertiary’ velocity field, which falls off like the inverse 
first power of the distance, and so on; the complete series for the remote 
influence of the obstacle is non-uniformly convergent, and the true behaviour 
which it represents is probably of a damped-wave type. ‘lhe remote in- 
fluence on a weakly-sheared oncoming flow of ‘ half-body’ shapes like the 
Rankine source body is also discussed. 

‘The method of this paper for calculating vorticity fields is an alternative 
to that developed by Hawthorne (1951, 1954) and Hawthorne & Martin 
(1955). Their integral for the streamwise component of vorticity is con- 
venient for numerical computation but was found unsuitable for analytical 
estimation near singularities and at infinity. ‘The present method is parti- 
cularly suitable for this purpose, and gives all components of the vorticity 
with equal ease. Detailed comparison between the results of the two 
methods for the sphere is given in §7. ‘This has been made possible by 
Professor Hawthorne’s kindness in putting at the author’s disposal data 
additional to those in the papers cited. 

One advantage of the present expressions for the secondary vorticity 
over those of Hawthorne (1954) accrues almost by chance. ‘The form of 
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these expressions, in the case of flow past an infinite cylinder, with the 
upstream velocity normal to its axis and increasing linearly with distance 
along the axis, is such that one can write down at once, by inspection, the 
secondary velocity field which the vorticity induces. ‘This is done in §4, 
and figure 1 shows the computed secondary flow for the case of a circular 
cylinder. This derivation of the secondary flow past cylinders appears 
somewhat as a sideline in this paper, but it may prove the most useful of the 
results which are obtained. 


2. THE VORTICITY FIELD IN WEAKLY-SHEARED FLOW PAST AN OBSTACLE 

Consider flow past an obstacle in which the velocity field far upstream 
is a parallel flow 

vz=V(y,3), =U, =0. (2) 
Suppose that the obstacle is near the axis y=2z=0, and that I’(y, s) varies 
only a little from its value ’(0,0)= U over a range of variation of y and 
large compared with the size of the obstacle. A particularly interesting 
case is that in which the gradient of V also varies little over this range, so that 
with suitable choice of axes we can take 

v,=U+-Ay, vy =v, =0 (3) 
far upstream, for some rate of shear 4. 

In all cases, the irrotatonal flow about the obstacle with uniform up- 
stream velocity U will be called the primary flow. The vortex lines, which 
far upstream lie in planes x= constant, are deformed by the primary flow into 
shapes which will be calculated, after which (§3) the resulting secondary 
vorticity field and the asymptotic form of the associated secondary velocity 
field far from the obstacle will be deduced. The work is displayed only 
for the simple case (3), but results for the more general case (2) are quoted 
at the end of § 3. 

Let the streamlines of the primary flow be represented by equations of 
the form 


V=H(X, Yor z= 2(X, Vo 2); (4) 
where Vo= lim (y), Z = lim (2). (5) 
t>— @ 


Thus expressions (4) are solutions of equations (1) for y and z, and can be 
supposed obtained by the ordinary methods of hydrodynamics. Then the 
solution for the additional variable ¢ is 
1 1 | 

where v(x, V9, 2p) signifies the value of the velocity component v, at the point 
x on the streamline given by (4). 

In case (3) the upstream vorticity field is given by 


w,=w,=0, w,=-~A, (7) 


so that the vortex lines are parallel to the z-axis. An element of vortex line 
stretching from (X, yo, 9) to (X, 9, 39 +539), Where X is large and negative, 
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is deformed by the primary flow, but always remains on a surface ¢ = constant 
and joins points on the streamlines specified by vp, 39 and Vo, 2+ 525. 
When its length is dr, the vorticity is changed from its initial value — 4 
to — Adr 5z,, and its direction lies along the new position of the element. 
Its components are therefore 


Ox oy Os 
(=) Vo (=). Yo (=), Yo ) 


where the suffixes indicate the variables to be kept constant in each 
differential coefficient. 


Now by (6), at ot ot 
= — dy + + dt, ) 
dt dVo dz, (9) 
so that Ox Ot Oz ot 
=-v,—. 0 
Yo ot Ox ) 


We have also =| (=) 
02 O39 ov Yo,20 03 
ov ot 


which, with a similar result for dz 029, gives by (8) 


ot ot ov ot Os 
=Av,—, w,=Av,— —A—, w,=Av,— -A 12 
Far downstream of the obstacle, we have 
x : 


where X(y, s) is Darwin’s ‘total drift,’ the ultimate displacement of a 
particle of fluid (relative to particles in its plane which are far from the path 
of the obstacle). We shall need Darwin’s result 


| | X(y,2) dydz=V,, (14) 


where I’, is the volume of fluid whose mass is_ the ‘hydrodynamic mass’ 
associated with the body’s motion. 

In the limiting conditions (13), the vorticity components (12) satisfy 
w, —> 0, —A, (15) 
so that the only new component far downstream is w,, which Hawthorne 
calls the ‘secondary trailing vorticity’. 

If the primary flow has a velocity potential U(x+¢), so that 4 is the 
‘disturbance potential’, there is a useful approximation to the solutions (4) 
and (6) of equations (1), which is valid on any streamline that remains far 
from the obstacle (so that (v,—U)/U=éd/éx<1 all along it). Such a 
streamline takes the form y = yo, ¥ == 2, toa first approximation, and toa 
second approximation 


wo, >A 


(16) 
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while equation (6) becomes 


Note that the integrals in (16) and (17) need not be taken along the stream- 


lines ; to the approximation involved, vy and s may be taken constant therein, 
and equal either to yy and sy or to their values at the point under consideration. 
Substituting from (16) and (17) in (12), we obtain to the same approxi- 


mation 


w,+A = —A | 


;quations (18) are valid far from the body, except downstream of it on 
streamlines which have not remained far from the body. On these, in fact, 


(15) holds. 


3. ‘THE SECONDARY FLOW IN WEAKLY-SHEARED FLOW PAST AN OBSTACLE 
The secondary flow associated with the secondary vorticity field is 


calculated in three parts: 


(i) a part, which is the gradient of a potential ¢, vanishing at infinity, 
4 and whose normal velocity component on the body surface cancels 


that due to the shear flow v7, = Ay, 7,=7,=0 which results from the 


undisturbed vorticity distribution (0, 0,— 4); 


(iii 


(ii) the Biot-Savart velocity field, say v,, of the vorticity change 
Ww, =(w,, w,, w, +41); note that the Biot-Savart integral would not 
converge if the undisturbed vorticity field were not subtracted out ; 
a further irrotational part, which is the gradient of a potential ¢, 
vanishing at infinity, and whose normal velocity component on the 
body surface cancels that due to {/:); this may be regarded as the 


velocity field of the image vorticity associated with w,, while part (1) 
is that of the image vorticity associated with the undisturbed vorticity 


field. 


Now, by general properties of harmonic functions, the irrotational 
velocity fields (i) and (iii) fall off like the inverse cube of the distance from 


the body, because the required normal velocity on the surface has a zero 


integral over the surface in each case (so that in‘ Green’s equivalent layer’ 
the total source strength is zero). However, the Biot—Savart field (ii) falls 


off more slowly, and its asymptotic form will now be derived. 
First, we see by inspection that a velocity field corresponding to the 


asymptotic form (18) of the vorticity change w, is 


f= dd 
v,=—A dx v,= Ad 
; (Its curl is (18), and its divergence vanishes). Since the disturbance poten- 


tial ¢ for flow past any finite body is asymptotically that of a doublet, this 


velocity field (19) falls off like the inverse square of the distance from the body, 


more slowly than that of (i) and (iti) above. 


(17) 
“2 a 2 
ws-A| — dx. (18) 
Os 
| 
= 
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Next, if w, is the difference of the vorticity change w, from its asymptotic 
form (18), and v, is its associated (Biot-Savart) velocity field, v, can be 
expressed asymptotically as the velocity field of a single vortex element of 
strength equal to the total strength of the field w,, thus 


where (x? + "+ 3?) and the integral is over the whole region occupied 
by the fluid. This step is permissible only after the inverse-cube part of w, 
has been subtracted out, as the integral would otherwise be only semi- 
convergent. * 

Now |, dz does not vanish, because of the trailing vorticity. If the 


integral is taken over a large rectangular box, we have 


| w,dr=| {div(ven,), div(ye,), div(se,)}d7 


= | (xv, dS 


= | | (x,y, 2)A 

=(0,0,-—AV,), (21) 
where (15) and Darwin’s result (14) have been used. Note that only the 
rear face of the box contributes to | , because the part of w, which falls off 


like 1 r* has been taken out. By (20) and (21), 
~ 0) as r— 0. (22) 


Now, for a body so placed that the irrotational flow exerts no yawing 
moment on it, the asymptotic behaviour of the disturbance potential (‘Taylor 
1928) is 


(23) 


where I’, is the volume of the body. Hence, for such a body, combining 
the results (19) and (22), we obtain 
A(V,,+2V,) 
vy ~ 0). (24) 
47 

More generally, é contains terms in y 7? and = r? proportional to the moments 
on the body, about the s-axis and y-axis respectively, and these make an 
additional contribution, through (19), to (24). It may be noted, however, 
that they do not affect the value of 7,,,0n y= = 0, which alone influences the 
first-order displacement of the stagnation streamline. 


* As a matter of interest, if this argument is incorrectly applied to the inverse- 
cube part of @, it gives an asymptotic contribution to the velocity field from this 
part only two-thirds as much as the true contribution (19). 
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It has been seen that a primary flow with the disturbance velocities 
falling off like ry * produces secondary flow velocities falling off like r-?. It 
may be questioned whether this conclusion does not undermine itself, since 
the assumption that the disturbance velocities are asymptotically given by a 
potential (23) is used to derive (24). However, the process by which (24) was 
derived from (23) was a linear one, so that the existence of a further con- 
tribution (24) to the disturbance velocities at infinity merely miakes an 
additional (call it ‘tertiary’) contribution to the vorticity field and the 
resulting velocity field. Calculations (along lines similar to those above) 
show that the asymptotic form of the tertiary velocity field is of order r 1, 
being in fact 

A*(V,,+2V,)(x? + y?) 


Srl 73 
Ty Us A*(V, +2V,). (3s? y*)(1+x r) (25) 


at all points far from the body except those immediately downstream of it. 

It is obvious that the series for the disturbance velocity whose first 
three terms are given by the primary flow L'V¢, the secondary tlow (24) and 
the tertiary flow (25), is of little value at distances from the body comparable 
with U A. In this region, where the disturbances are very small, a truer 
picture may be obtained by solutions of the linearized form 


ope curlv= — Ax (26) 
Ox Ox 


of Helmholtz’s equation for the vorticity. Solutions of (26) which tend to 
zero at large distances are usually wavy, for example 

V~Vor “exp[i(ax+ by U], (27) 
where a,b,c are numbers satisfying Evidently (27) 
expanded in powers of 4 U would give a series of the kind described above, 
but its convergence would be non-uniform for r large, and the true order of 
magnitude of the whole is no greater than that of the first term. 

The above gives a mathematical reason for not feeling alarmed at the 
orders of magnitude found. Physically, one would not expect the approxi- 
mation to be good at distances from the body comparable with that at which 
the oncoming flow differs from the primary flow by an amount equal to itself. 

We consider next how the results are to be modified if a body experiences 
drag, so that the irrotational flow outside the wake has a source-like behaviour 
at infinity (corresponding to the reduced mass flow in the wake). ‘The 
same problem arises for ‘half-bodies’ extending to infinity downstream 
such as pitot tubes—which exhibit the same source-flow behaviour. In 
either case the disturbance potential satisfies 


dbw—-— asr—> (28) 


where, for a finite body with a wake, m is the drag divided by pl *, while, for 
irrotational flow about a half-body, m is the cross-sectional area of the 


4 
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semi-infinite cylindrical portion. By (19), the secondary flow then is 
asymptotically 


Amx x Am 


far from the body except immediately downstream of it. ‘The contribution 
of the Biot—Savart field of w, falls off like »* and so does not have to be 
included in (29). As with (24), (29) is probably only a term in a series 
whose early terms are a poor approximation for values of r comparable with 
U/A. 

Finally, we give without proof some results for the general case of a 
parallel flow given by equations (2) far upstream. ‘The secondary vorti- 

where IV’, stands for l’(yo,3,) and the notation for Jacobians has been used. 
The secondary trailing vorticity is 


on 
The asymptotic form of the vorticity change w,, as y+ % except immediately 
behind the body (corresponding to (18), is 


od 
a( v) 


The asymptotic velocity field v, corresponding to (32) is less simple than 
one would like, because the expression analogous to (19), namely 
(- <4), (33) 
Oy) dy J dy dz 

although it has a curl equal to (32), is not solenoidal unless 
dy? +071 és?=0. In general, then, v, is asymptotically (33) minus 
grads, where V2;=4V2V. 

As in the special case, the secondary flow ts asymptotically v, + v,, where 
v, is the velocity field of a single vortex element of strength equal to the total 
strength of the field w,. We can show that 


Sox dvds) (- 5,0. 


> 

- 
= Od 
— dx, V 
~ 
1 
— 

| 3X 5,0), (34) 
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where each integral can be interpreted as the product of I’, (see (14)) and 
a mean shear over a region comparable with the body size. 


4, "THEORY FOR TWO-DIMENSIONAL PRIMARY FLOWS 

The theory for a uniformly sheared oncoming flow (3) will now be 
specialized further to the case of a two-dimensional primary flow in the 
(x,z)-plane about a cylinder with generators in the y-direction. (Note that 
the assumed oncoming vorticity (7) is entirely in the z-direction, but that 
the theory can be trivially extended to cases where there is in addition some 
uniform oncoming vorticity in the y-direction, since vortex elements in that 
direction are not stretched at all by the primary flow.) ‘The case in question 
is remarkable in that the exact secondary flow can be calculated. 

‘The streamlines of the primary flow may be written 


%=2(x, (35) 


and the time ¢ (see (6)) is a function of x and 2, alon«, giving, by (12), vorti- 
city components in the form 


él ot 0s 
Wy =. w, =9, w, = Av, (36) 
~0 


Some simplification is possible for these flows, however, if the time #, once 
calculated, is expressed as a function of x and zs, rather than of x and 2». 
Since ('s,, by (5), is the stream function of the two-dimensional flow, 


(37) 
Os tig dz 
Using these results to transform (36) for the case when f¢ is expressed as 
t(x, =), we obtain 
Ot _ Ot 
w,=AU— w,= 38) 
Os ’ y ( J 

Thus the vorticity is simply proportional to grad ¢ in magnitude, but 
is perpendicular to it in direction (as it must be because f= constant along 
vortex lines). "This result could also have been obtained by physical argu- 
ment: the rate of stretching of the vortex elements on which ¢= constant 
must be inversely proportional to the distance between successive lines 
t=constant by conservation of mass. 

‘The exact secondary flow can be deduced from (38) by inspection, in 
three parts (i), (ii) and (iii) as in § 3. The potential 4, whose gradient cancels 
out the normal component of the shear flow (Ay, 0, 0) on the body is 

= Ayd, (39) 


where A(x,3) is the ‘disturbance potential’ as in $§ 2 and 3. The Biot- 
Savart field of the vorticity change 


t 
(40) 
Ox 


ist v,=0, v, = A(x— Ut), v,=0, (+1) 


| 
‘ 
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because the curl of (41) is (40) and its divergence vanishes. Finally, no 
extra irrotational portion (iii) is needed because the velocity field (41) is 
wholly tangential to the body. 

Combining the secondary flows (39) and (41) with the primary flow and 
the oncoming shear flow, we see that the velocity components in any plane 
y=constant are simply 


v,=(U+Ay)(1+—}, v,=(U+Ayv)—, 42 
(42) 
namely, the ordinary potential flow with the upstream velocity appropriate 
to that plane. Superimposed on these motions, however, is a flow parallel 
to the generators given by 


A(x+o— (43) 


Figure 1. Primary and secondary flow about a circular cylinder of radius a, with 
axis the y-axis, when the upstream velocity is (U+ Ay, 0, 0). 
Streamlines of the primary flow (note that the velocity components in the x and 
z directions follow these streamlines even when the secondary flow is included). 
Contours of constant 7, Aa (note that the negative values, which predominate, 
denote secondary flow in the direction of decreasing primary flow velocity). 
(Flow is from left to right and only the upper half is shown.) 


As an example, figure 1 shows the distribution of v,/Aa for sheared flow 
about a circular cylinder of radius a, computed from Sir Charles Darwin’s 
evaluation of ¢ for a circular cylinder. Remembering that the y-axis is 
along the axis of the cylinder in the direction of increasing velocity, we see 
that the greatest secondary flow velocities are in the direction of decreasing 
velocity—physically, because this is the direction of decreasing stagnation- 
point pressure. Like ¢ itself, the secondary flow velocity has theoretically 
a logarithmic infinity on the surface of the body (and on the central stream- 
line behind the body); but viscous resistance must place a practical limit 
on the velocities which can be achieved. 
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5. THEORY FOR AXISYMMETRICAL PRIMARY FLOWS 


Some simplifications, analogous to those of §4, are possible when the 
primary flow is axisymmetric, though unfortunately they do not extend so 
far as a simple exact solution for the secondary flow. 

If x, p, A are cylindrical polar coordinates relative to the axis of symmetry, 
the streamlines of the primary flow may be written 


p=p(x, py), A=constant. (44) 


In ordinary hydrodynamics such streamlines are obtained in the form 
w%=constant, where y is Stokes’s stream function; note that the relation 
between us and p,=lim (p) is 


b=4Up,?. 


The time l 
t=t(x, po) = + — 46 
Py) | v(x, Po) ( ) 
is independent of A. 
If p and A are related to the coordinates x, v and z used in earlier sections 
by the equations y=pcosA, s=psinA, then the upstream flow (3) and its 
associated vorticity field may be written 


=—Asind, w,=—Acosd. (47) 


v,=U+Apcosa, v,=7,=0, w,=0, w 4 


‘The component w,, the ‘ring vorticity’, varies in a particularly simple 
manner due to the stretching of vortex elements by the primary flow. An 
element of ring vorticity subtends always the same angle 6A at the axis, and 
so its length varies in direct proportion to the radius p. At any point, 


therefore, 
) U 
> (—A cosA)£ =(—AcosdA)p J): (48) 
Po 


where (45) has been used to throw the result into an alternative form in- 
volving Stokes’s stream function. 

The radial and axial components of vorticity, w, and w,, are stretched 
in a manner depending more closely on the details of the flow in planes 
A=constant, and arguments precisely similar to those leading to (12) give 
for these components the results 


Ct Ct = Op 
w,=(dAsinA)v,— , w,=(A sind)(e (49) 
Accordingly, the secondary trailing vorticity is 
lim (w,)=(A sinA)X(p). (50) 


Since, as (46) shows, ¢ is necessarily calculated in the first place as a 
function 1(., p,)), the formulas (49) in terms of that function will be found 
most convenient for practical application in §7. It is interesting and 
instructive, however, by analogy with §4, to see what they become in terms 
of a t which, once calculated, is expressed as a function of x and p, rather 


: 
— 
. 
3 
5; 
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than of vandp,. From (45) and the expressions for the velocity components 
in terms of Stokes’s stream function, 


Ox po) Op po] 


Using these results to transform (49) for the case when ¢ is expressed as a 
function t(x,p), we obtain 
(52) 
p Op p Ox 

‘Thus the vorticity resultant in each meridian plane, while perpendicular 
to grad ¢ as in $4, is proportional in magnitude to (pp) grad ¢. This, like 
the corresponding two-dimensional result, can alternatively be derived 
from conservation of mass. ‘The mass of fluid in an annular region bounded 
by two neighbouring axisymmetrical streamtubes and by two neighbouring 
surfaces ¢=constant is proportional to the product of the radius p, the 
distance between the two surfaces t = constant, and the length of an element 
in which the part of a surface t= constant between the two streamtubes is 
cut bya meridian plane. ,For this product to remain constant along a stream- 
line, the length of such an element must vary as (grad t) p, whence the 
required result follows. 

The author could not obtain by inspection the velocity field corres- 
ponding to the values of w,, w, and w, given by (48) and (52); a field whose 
curl is as required is fairly easily found, but it is not solenoidal, nor does it 
satisfy the boundary conditions! ‘Therefore, the direct use of the Biot- 
Savart formula for calculating the secondary velocity field seems to be 
necessary. 

Finally, it may be noted that in some problems it is convenient to use the 
spherical polar coordinates r, 4, where 


wy = —(Asind)U 


rcos@= x, rsin@=p, (53 
instead of x and p. ‘The ‘longitude’ A is retained. | By considering the 
stretching of a vortex element joining streamlines specified by (py, A) and 
(py + dp9, A) we obtain (exactly as (8) was obtained) the equations 

w,=(—Asind) (=) A sind) ( r (54) 
Opa) \ ts 
where equation (47) for the initial vorticity of the element has been used. 
If the time is expressed as t(py, 7), then 


cr it Ot cu taf Or 
Opo/ t.i ot or OPpy po/ 
as in (10) and (11), giving 


w,=(AsinA)zv, —, w,=(Asina =) (56) 
"Cpy \ ” Opy 
where 7, .\ are kept constant in the differentiations. If, alternatively, the 
time is expressed as ¢ (p,,@), then 


au Cre Ct or er c.f 
& = {— + r— , (57) 
t.2 ot ret CPo/ 1.3 Opo/ 6,3 Ca\ 


= 
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=(A sin dA)z w,=(AsindA) (: (58) 
where 6, A are kept constant in the differentiations. If, lastly, the time ts 
expressed as ¢(7,#), then the sentence following equation mes makes it clear 
—- = w,y= —(AsinaAyU = (59) 


6. EVALUATION OF t FOR FLOW PAST A SPHERE 
The drift function ¢ will now be studied in detail for flow past a sphere, 
principally so that the results can be substituted in the formulas of §§2, 3 
and 5, and the secondary flow evaluated. 
Spherical polar coordinates r, 4, A are used. From Stokes’s stream 
function for flow past a sphere, and equation (45), the streamlines are given 


by 


where a is the radius of the sphere. Two alternative equations for 
obtaining the time ¢ are useful in different regions, namely 


dr dr 
U (1 -os 
and dt - (62) 
v(t sin @ 


On any given streamline, by (1), 
t+rU—-0 as0--7, (63) 
since v+r -0. ‘The integral for ¢ derivable from (60), (61) and (63) is 
unfortunately hyperelliptic, as Darwin (1953) points out. However, 
valuable expansions of it can be obtained both for large and small values of p, a. 
For large p, a, the solution for r of equation (60) can be expanded in 

powers of a® as 
sin 4 S pi 


It tollows that (62) can be expanded as 


12 
55) 


sin? 
+ 

2r3 

3 5 712 

ode 878 167° 12872 
sin?4 


=— —, (1+ 5 —sin*+ — (05) 
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Integrating, 


~— | Si dé— sin‘é dé+ | sin 
8 pis Po 128 pj! 

(66) 


where the integrals have been made definite in such a way that (63) is satisfied. 
The radius of convergence of the series (65) is easily determined from the 
position of the singularities of r/(1+ as a function of py cosec 
The series converges for py ~3!72 1%asin@. Hence (66) converges for all 
@ provided that pp, a > 322-13 = 1-375*. 
Darwin’s drift distance is inferred from (66) as 


3 (37 (32 315 a? 
6 9 12 
=0-442 0-914 -..., (67) 
Po Po Po 


the fractions in brackets being the values of the integrals in (66) for @=0. 
Note that (66) has the properties (obviously true on any streamline 
from the symmetry of the flow) 


—t(47) =t(47)— (7-8), 
J 


For small pg a, each streamline (60) is divided into two parts, one on 
which @ is near 0 or z (that is, sin @ is small) and one on which (r—a) a is 
small. When @ is near 7, we use the approximation 


—sec@ = 1+ }sin?6 (69) 
in (61), together with (60) and (63), to obtain 


a~ 


\ (7? +ar+a?) 


1 Por 1 ay 3 


We suppose that this approximation is valid down to r= 7, (say), after which 
@ departs too greatly from 7 but (r—a) @is small. For r<ry, we use (62) 
with the approximation 


2 4+ py” 


@ 9 


* Some improvement in convergence arises if the series (66) is rearranged as a 
series in powers of a/r (with coefficients functions of 6), using (64). This has not 
been done, however, because of the desirability of knowing t as a function of po for 
the application of the results of §5. The results it gives for X(p) are identical with 
those obtained by the method used below. 


Ut=p,coté+ 
. 
1 
; 
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(in which the error is of the order p}) to obtain 


Ut= 5X (po) + 


log cot} 4+ cot cosec (72) 
ai} : Ya 


where the arbitrary constant is determined by putting 4 = 57 and using (68). 
Comparing (70) and (72) for values of r and 4 for which both approximations 
are valid (and also (60) holds), we obtain after some calculation 


X(po) = 3 a(1 (2 (- 3 1) (73) 


As Darwin points out, \(p,) is logarithmically infinite at pyp=0. This is 
because there is an infinite time-delay as fluid approaches any stagnation 
point where the tangent to the body surface is smoothly-turning. 

In figure 2, the function pX(p) a? is plotted against p a, the expression 
(73) being used for p a <0-4, and expression (67) for p @> 1-5 (some allow- 
ance for the later terms in the series being made up in this case by assuming 
that the coefficients remain in approximate geometrical progression). ‘T’o 
interpolate between them we use Darwin’s result (14) to give 


| 72 | X(y,3) dvds (; 7a ) 


(74) 


on 
since the hydrodynamic mass of a sphere is half the displaced mass of fluid. 
Practically no choice is possible if a smooth interpolation between the two 
crosses in figure 2 is to be made which satisfies (74), and a ‘ French curve’ 
has been selected and used for this purpose which makes the integral from 
0-4 to 1-5 (calculated by Simpson’s rule) make up exactly for the amount by 
which the sum of the (analytically evaluated) integrals from 0 to 0-4 and 1-5 to 
© falls short of }. The values of X(p) a itself as a function of p a have been 
read off figure 2 and are presented in figure 3 and table 1. 

Once X(p) has been obtained, it is a simple matter to compute the function 
t throughout the field. For #87, (69) is a good approximation, and 
accordingly the formula (70) has been used for ¢. For $7<0<27, two 
formulas have been used: when r,a< 1-5, (71) is a good approximation, and 
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Figure 2. The method of interpolation used to derive Darwin’s ‘total drift’ 
function X(p) for irrotational flow past a sphere. 
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Figure 3. Darwin’s ‘total drift” function X(p) for irrotational flow past a sphere. 
Since the scales used for p and X(p) are the same, the figure gives the actu 
shape into which planes of fluid, initially at right angles to the stream, would 


ultimately be distorted in such a flow. 
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so formula (72) has been used for t, while when ra -1-75, we have 
p, cosec # 1-58, which is well inside the radius of convergence of (65), so 
that (66) has been used. he interpolation across the small gap is straight- 
forward. For @- }z, equation (68) has been used to deduce ¢ from its 
values when # -\>. ‘The results are exhibited in figure +, where the surfaces 
t=constant are shown for values of Uta proceeding by intervals of 0-5 
from —2:5 to +3. The streamlines p, a=0-125, 0-25, 0-375, 0-5, 0-75, 1-0, 
1:5 and 2-5 are also shown. 


| 
| 


| 


Figure 4+. The irrotational flow past a sphere : streamlines and surfaces ¢— constant. 
The latter are the shapes into which planes of fluid initially at right angles to 
the stream would be distorted as they passed over the sphere (if the flow 
were Irrotational) 


Figure + has two uses. First, it exhibits a steady fluid motion in far 
more detail than has ever been given before, showing precisely what happens 
to every particle of fluid and when (Darwin’s figure 1, representing flow past 
a circular cylinder, is insufficiently detailed for this purpose.) Secondly, 
in the problem of shear flow past a sphere, it exhibits precisely how vortex 
elements initially in the p-direction are stretched and change their direction 
during passage of the sphere. An impression of the values of the vorticity 
resultant in a meridian plane can be derived either by observing this 
stretching, or by means of the rule derived in $5 that along any streamline 
the product of this vorticity with the distance from the axis varies inversely 
as the distance between successive surfaces t=constant. Numerical 
values of the vorticity components are obtained in § 7. 


| | | | | 
3 | 
} | 
| 
| | 
| | 
| } | 
| | 
| | 
| | | | 
| | 
| 
—— 
} | 
| | | | | | 
| | 
| 
| 
(ff 
: 
\c 
f 
\ 


Drift 49 


An alternative mode of display of the function ¢ is adopted in figure 5, 
which shows the variation of Uta along all but one of the streamlines in 
figure +, by plotting Uta asa function of «=rcosé for each value of p,/a. 
A third manner of representation of some of the results (figure 6) will be 
discussed in § 7. 
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Figure 5. Variation of the drift function t, for flow past a sphere of radius a, along 
different streamlines, which are identified on the figure by the values of po a 
noted on each curve. Here, po is the distance of a streamline from the axis 


far upstream. 


7. THE VORTICITY FIELD IN WEAKLY-SHEARED FLOW PAST A SPHERE 

The drift function ¢ having been evaluated for irrotational flow past a 
sphere, equations (56) and (58) show that only its derivative with respect to 
py (keeping either 7 or @ constant) is needed to deduce both w, and w, in the 
problem of weakly-sheared flow past a sphere. 

In the region 6 > 27, where ¢ has been expressed as a function of py and 


“6 
r as in (70), equations (56) are applicable; in applying them, we note that 
by (60) 
rtané@ (75) 
r({—) = ‘ 75 
Opo r Po 
In the region }7<@<37, where t has been expressed as two different 


functions of p, and 4 in different ranges of py, the detailed computations of t 
were done for 6=90°(10°)150°, with a suitable selection of values of py in 
each case ; a small gap is present (1-5 <r @<1-75) where neither formula is 
reliable. However, as figure 6 shows, a smooth curve can easily be drawn 
through the points which are available for each 6. From these (@t/dpy), 
can be inferred by measuring the slope of the curve ; for the larger and smaller 
values of py a, however, analytical differentiation is preferable. The 
resulting tables of values can then be improved by the technique of 
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Figure 6. Variation of the drift function ¢ with py for constant values of the spherical 
polar coordinate 
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‘graduation’, or smoothing of data. When this has been done, w, and w, 
are inferred from (58), in which, by (60), 
pd 3 
(5-) (76) 
Opo/s + 

‘The same technique applied to the values of X(p) itself gives the important 
function X’(p), which by (50) determines the secondary trailing vorticity 
distribution. Since X’(p) - as p 0), it is convenient to exhibit the 
product —pX"(p) in Table 2. Certain of the values can be compared with 


pa ( 5 0-1 0-15 0-2 25 0-3 0-35 
pX'(p) a | 1-333 | 1-325 | 1-307 | 1-283 | 1-253 | 1-214 | 1-166 | 1-11 


pia 04+ |06 [0-7 [08 [1-0 4-1 
~ pX'(p)/a | 1-043 | 0-897 | 0-749 | 0-621 | 0-523 | 0-450 | 0-394 | 0-345 
pla 1-3 1-4 1-75 | 2-0 
| pX(p)/a | 0-298 | 0-249 | 0-198 | 0-151 | 0-079 | 0-049 | 0-030 | 0-019 


Table 2. 


results for the secondary trailing vorticity given by Hawthorne & Martin 
(1955), who obtain values of (w,),_,... d sind equal to 12-976, 4-792, 1-852 
and 0-404 for py a@=0-1, 0-25, 0-5 and 1-0 respectively. ‘These values, 
obtained by numerical integration along streamlines, differ from those to be 
obtained from ‘Table 2 (by dividing by pa) by —0-09, —0-06, +0-06 and 
+0-01 respectively. Errors in both methods probably contribute about 
equally to these discrepancies, although in the case p a=0-1 the value in 
Table 2, based on expression (73) for small p a, can perhaps be preferred*. 

The general distribution of vorticity is shown in ‘Table 3, where the 
variation of w,, and w,, along three particular streamlines py a=}, } and 1 is 
exhibited by giving their values at 10° intervals of the polar angle 6. ‘These 
streamlines were selected partly because they give a good general idea of the 
variation of w,, and w,, over the field of flow, and partly because they give a 
convenient basis (as will be seen in a later paper) for the evaluation of the 
Biot—Savart field of w,, by integration first with respect to A, secondly with 
respect to @ along a streamline py = constant, and finally with respect to py. 
The vorticity change w, is tabulated as being the quantity directly respon- 
sible for the secondary flow. ‘lhe absolute vorticity components w, and w, 
can, however, be obtained by adding to w,, and w,, the values of the undis- 
turbed vorticity components (w,= — A sinAsin@, w,= — A sind cos@) which 
are tabulated on the left. 

Professor Hawthorne has kindly supplied the author with the results of 
computations, additional to those given in Hawthorne & Martin (1955), 
from which values of w,, and w,, can be deduced. Comparisons with the 
values to Table 3 at a sample of points show agreement to within the same 
order of accuracy as was noted above for the secondary trailing vorticity. 

*_— X'(p) ‘a is given as a series of powers of (p/a)* of which the first two terms 
tor p/a=0-1 are 1°333—0-026 ; the next would be expected to be less than 0-001. 
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For values of py @ which are either large or small, the direct use of 
equations (56) and (58), with the approximation (66) or (70) or (72) for ¢, is 
recommended in preference to the numerical procedures described above. 
Here, to save space, we give only the first approximation to w,, and wy4, in 
each case. For large r a, except immediately downstream of the sphere, 


w,, ~ AsinAsin#é(-} , wy ~ —sAsinAcosd(-} . (77) 
r r 


For 6 near =, 
( 3\12 92 3\ 3/2 
wy, ~ AsinAsind (1-5) 


a+2r 


~ 


8) 


( —1/2) 
wy ~A sindeos@{1-(1~ 5) 
J 


For 6 near 0, we must add on to the values (78) the considerably larger 
values 


w4,~A sindeos (1 )X | 
(79) 
wip ~ —AsinAsin a(1 | 
For small (r—a)/a, 


Finally, we note the simple exact form of the third component w, of the 
secondary vorticity, which by (48) and (60) is 


3\—-1/2 3\—-12 


It is interesting that this component shows no dependence on @. 
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On the thrust due to an air jet flowing from a wing 
placed in a wind tunnel 


By L. C. WOODS 


University of Sydney, Australia 


(Received 20 December 1955) 


SUMMARY 
Consider a wing-jet combination placed in a wind tunnel; the 
measured thrust on the wing due to a high speed jet emerging from 
it will require to be corrected to give the infinite stream value. 
This paper provides a theory of these wind tunnel corrections, and 
incidently establishes that the ideal thrust is almost independent 
of the jet exit angle. 


1.. INTRODUCTION 

The figure shows a section through a model of a wing and jet engine 
combination, the outer and inner surfaces of which are S, and S, respectively. 
This model is in a cylindrical wind tunnel, S,, of arbitrary cross sectional 
shape. Air is drawn into the wing at BBL’, the surface S, separating the 
external stream from that drawn into the engine. It is assumed that energy 
and mass is added to this stream at some section EE” within the wing, con- 
verting it into a high speed jet which emerges at some angle to the main 


Wing-jet combination in wind tunnel. 


stream at the trailing edge DD’. In our idealized mathematical model ot 
the flow this jet is separated from a dead-air region (shown shaded in the 
figure) by a vortex sheet S,, and a further vortex sheet S, separates the dead- 
air region from the external flow. The velocity in the dead-air region. 
which is introduced to represent the blockage effect of the wing drag, is 
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assumed to be zero. ‘Thus the cross sectional area of this region is equal to 
the displacement area produced by the action of viscosity over the wing 
surface and internal ducts*. 

Let the plane 1, L’ , cut the tunnel at right angles upstream at infinity, 
and let the curves of intersection of this plane with the surfaces S, and S, 
enclose areas 4 and a respectively. Similarly, let the surfaces S, and S, 
have cross sectional areas (b+ /,) and 4, downstream at infinity, so that 5 is 
the displacement area mentioned above. 

So far we have assumed in our idealization of the flow that (1) there is no 
mixing of the jet and main stream fluids, and that (2) the vorticity is con- 
fined to vortex sheets. With these assumptions it is natural to assume 
further that (3) the flow is isentropic in both the jet and the main stream. 
‘Thus our model is only a crude first approximation to the real flow, but it 
does yield what might be termed the “ideal thrust”, and a comparison of 
this with the actual thrust achieved will give a valuable estimate of the loss 
of jet efficiency due to viscous and turbulent action. 

Values of the velocity U’, pressure p, and density p at infinity, (a) upstream, 
(6) downstream in the external flow, and (c) downstream in the jet, are dis- 
tinguished by the subscripts ,, , and , respectively; the subscript , is 
used to denote values in the jet at the trailing edge. The boundary con- 
ditions are (1) no flow across tunnel and wing surfaces, (2) continuity of 
pressure across the dead-air region and the vortex sheets S, and S,, (3) 
continuity of pressure and velocity across the surface S,, and (4) known 
sources of mass and energy within the wing. By (2) the pressure down 
stream at infinity is p, right across the tunnel, i.e. p) =p,,- 


2. EXACT MATHEMATICAL THEORY 
Let n be a unit vector along the outward normal to the (external) region 
bounded by the surfaces S,, S,, S,and the planes L_, L’ M, M3, 
and n’ a similar outward vector for the (internal) region ©, bounded by S,, 
S, and the planes L_,L’,, M, Let i be a unit vector directed 
downstream along the tunnel axis. ‘lhe upstream thrust 7’ acting on the 
wing may then be written 


T=- pn.i dS — | pn’.i aS (1) 
J S3+S; 


*Alternatively (see Allen & Vincenti 1944) this cross sectional area can be chose? 
so that it gives the same pressure drop downstream at infinity as that caused bv the 
drag in the actual non-isentropic flow downstream of the model. ‘The displacement 
thickness method, which is equivalent to assuming isentropic flow, is preferable in 
the author’s opinion, since the non-isentropic flow of the wake will not fill the tunnel 
for a considerable distance downstream of the model (many tunnel widths), whereas 
the lower pressure drop of the isentropic flow will be almost completely achieved 
within a downstream distance of less than one tunnel width. The pressure gradient 
and increased velocity at the model due to the wake blockage will be produced pre- 
dominantly by this latter pressure drop. 
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The right-hand side of this equation contains no drag force, the dead-air 
region representing only the displacement effect of the drag. ‘lo represent 
the drag force fully we must also remove momentum from the stream. 

If q is the velocity vector and S is a surface in the fluid with a normal n 
the momentum equation can be written 


| pn dS =—| (n.q)pq dS. 


We apply this equation to the regions %, and ¥,, and take the component in 
the i-direction to obtain 


| pn.i dS —(A—a)p_, +(A—b—hy)px 
- 
| dS ~ap_, +hops = ap_.U®, —hypol 
On adding these ee and using (1), we find that 
hg(py U5 — — A(p_ + (2) 


In this equation the upstream values and the areas 4 and 4 (see later) are 
known quantities, and the six quantities p,, p,, U,, py, Uy, My have to be 
calculated. ‘The six equations necessary for this are 


(3) 

ap_x l =hy pol (4) 

y-1 pow y-1 Px 

—1p, po 


in which Q’ is the known strength of the source within the wing and the 
numbers y and y’ are the ratios of the specific heats in the main stream and 
jet respectively. Equations (7) and (8) are adiabatic equations of state, 
while the first and second pairs follow from i conservation of mass and 
energy, respectively, for each of the regions X, and = 

Equations (2) to (8) contain the exact solution ie our problem, but for 
practical purposes the following approximate theory is sufficient. 
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3. APPROXIMATE THEORY 
For wind tunnels of cross sectional dimensions large relative to the 
model dimensions, we can set 
U,=U_,(1+8), (9) 
where 6 is a small first order term, and disregard terms of the second order 
in 6. Of course, in the limit as the tunnel cross section tends to infinity, 
6+ 0, and the approximate theory will yield exact results for an infinite 
stream. From (5), (7) and (9) it is not difficult to calculate the expansions 
p,=p_,{1— + [(2 —y)M?— 1] + 0(63) }, (10) 
and 
U2 ,3(1 + + 0(8%), (11) 
where VV is the Mach number upstream at infinity, and B=(1—M?)!?. 
With these approximations, equation (3) gives 6?6=(b+h)—a)/A, to the 
first order in small quantities, and 


A 2p! 
to the second order. Equations (10) to (12) enable us to replace (2) by 


where c is the planform area of the wing. 
The areas a, 6, c and h, can be eliminated from (13) by introducing the 
following coefficients : 


the thrust coefficient 


T 
Cr= (14) 
the momentum coefficient of the jet at infinity, 
hopo U5 


(15) 


~J0 
Lep_ 


and a similar coefficient C,, referring to conditions at the jet exit ; 
the mass coefficient of the air intake, 
ap_U_. a 
the mass coefficient of the source within the wing, 
Q cp_,U ( ) 
and the drag coefficient, 


4 
(18) 


Cp= 


where D is the drag, which is equal to bp_,U?,, the rate at which 
momentum is destroyed. 
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The momentum coefhcients C,), and C,, can be related by using the 
equations (6), (8), and hy py Uy =h, p, Uj, the last of which expresses the fact 
that mass is conserved in the jet. Thus we find that C,,=C,,(U, U); and 
if the pressure change along the jet is assumed to be small, 


} 0 1 


from the definition of C,, and the fact that p, =p... However, for simplicity 
we will express our results in terms of Cy). The drag coefficient with the 
jet operating may differ appreciably from that for the wing alone. ‘The jet 
will induce high speeds in the external flow near the trailing edge and will 
probably prevent any flow separation, but, on the other hand an additional 
boundary layer will emerge from the jet duct. It will not be easy to deter- 
mine experimentally the appropriate value of C’, for (18), but as a first 
approximation the value of (’',, for the wing alone could be taken. 
From equations (4), (15), (16) and (17), it is found that 


h,, = 2c(€ Q al Jo (20) 
where o=py p_, should be near enough to unity in practice. Equations 
(13) to (20) now give the basic result of this paper, namely 


In the limit 4 - ~, which corresponds to the case of an infinite stream, 
(22) 


Equations (21) and (22) apparently show that (’, is independent of the 
jet exit angle and the geometry of the wing but this conclusion is not really 
true. ‘The momentum coefhcient that must be adopted for practical 
purposes is the trailing edge value (',,, and if this value is maintained 
constant, equation (19) shows that (’,, will depend on the pressure p, at the 
trailing edge. In turn, p, will depend on the wing shape, and especially 
on the jet exit angle, 7, say. ‘The determination of the precise dependence 
of p, on + would demand a much more detailed and accurate model of the 
flow than the one we have taken. However, it seems clear from (19) that 
(, will not vary so rapidly with 7 as the function (1 — cos 7)—the law of 
dependence our engineering intuition at first suggests. ‘The independence 
of ©, and + does follow it (, is assumed to be constant along the jet. This 
was first noted by Davidson (1955) for the two-dimensional case. The 
practical importance of (22) is that, by deflecting jets downwards, lift might 
be obtained without much corresponding loss in forward thrust, but of 


course the effects of turbulence and viscosity have vet to be considered. 
Further, the induced drag could well exceed that produced when lift is 
obtained in the conventional way by wing incidence. High lift at low 
speeds rather than ‘inexpensive’ lift would seem to be the real advantage 
of jet deflection as a lifting device. 
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‘Three special cases of (22) are of some interest. (1) When the jet is derived 
from a source within the wing alone (as in rockets), a=0 in (13), and (22) 
reduces to 

(2) When there is no jet, but air intakes are operating (as in suction on wings 
for boundary layer control), (22) gives 
Cr 2C 

so that there is a ‘sink-drag’ 2C, on the wing. ‘This is a well-known 
result. (3) A special case of (1) occurs when the jet has the same total head 
as the main stream. From (i5) and (17), we find that, in this case, 
(,,=2C/,, and so there is a ‘ source-thrust’ given by 


@ 


4+. WIND TUNNEL CORRECTIONS 

Except for the wake blockage term, the theory of wind tunnel inter- 
ference on the thrust is almost the same as the theory of the interference on 
the drag of wings. (For a summary of this latter theory, see Howarth 1953, 
p. 522.) 

Let p*, be the free stream values corresponding to U p_ 
such that the velocity at a given point on the wing is the same in the free 
stream as in the tunnel stream. Let €,, €,, be the solid and wake blockage 
factors respectively, that is 


‘he factor «, depends on the tunnel cross section and model volume, and 
its values for two- and three-dimensional tunnels have been published 
(Howarth 1953). On the other hand, «,, is independent of the tunnel and 
model shape, being equal to one half of the value of 6 defined in equation (9). 
(A semi-infinite disturbance like the wake of a wing always progluces half 
of its final downstream effect of the wing itself.) Thus, from (12), 


= 
Ut, =U (23) 
Subtracting (21) from (22) we have 
Ce 


where the stars denote free stream values, and k is the last term in (21). 
Thus, for fixed values of jet momentum and intake mass, 


CF = (Cy— 2) 1) 


which, by an equation for p*,, corresponding to (10), and by (23), we can 
write as 


h+h,—a 


= 

q 

\ 

/ | 

a 

Be 
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By (16), (18) and (20), it now follows that 


aC 
(CytCo)? ] 
482.4 | D aC ( 4) 


which shows how the wind tunnel value Cy, may be reduced to the free 
stream value C7. 
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Vibrational relaxation in oxygen and nitrogen 
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Palmer Physical Laboratory, Princeton University 
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SUMMARY 

A converging channel of area ratio 34 : 1 has been developed to 
produce strong shock waves in a tube. Shock waves of speeds 
M,=3-7-5 in oxygen and W,=5-—10 in nitrogen have now been 
studied with an interferometer, and values of the relaxation time 
+ for the approach to vibrational equilibrium behind the shocks 
have been measured. ‘The value of 7 (atmospheric) varies from 
54 usec at 800° K to 1-3 usec at 3000° K for oxygen and from 19 psec 
at 3500° K to 5 usec at 5500° K for nitrogen. For oxygen, the graph 
of logz against 7-13 is not quite the straight line predicted by the 
Landau—Teller theory. The density ratios across the shocks were 
measured and compared with values calculated by the Bethe—Teller 
method for variable specific heats. Agreement between the 
measured and calculated values is satisfactory. Experiments 
were also performed on oxygen—nitrogen mixtures to determine 
the effect of nitrogen on the approach to equilibrium of the oxygen. 
It was found that O, and N, collisions at approximately 2000° K are 
40°, as effective in transferring energy to the oxygen as O, and O, 
collisions. A device that detects a shock with a time lag of less 
than 1 psec, consisting of an evaporated gold film which changes its 
resistance when heated by the shock, was also developed. 


INTRODUCTION 

‘The purpose of this investigation was to determine the state variables 
for a strong shock wave whose Mach number, VV, =7, 4, (7, is the velocity 
of the incoming gas relative to the shock and a, is the velocity of sound 
in it), is so high that the approximation of constant specific heats is no 
longer applicable, and, where possible, to measure the relaxation times for 
the approach to equilibrium behind the shock.* Particular etfort was made 
to study oxygen and nitrogen in the range V7, =3-10, and the results are 
compared with certain theoretical work concerning relaxation effects, and 
also with the relatively scarce experimental data for O, and N, that already 
exists. 

* A fuller description of the work reported in this paper can be found in the 
author’s doctoral dissertation, on file at the Princeton University library. 


61 
: 


62 Vernon Blackman 


Owing to variation of the specific heats of oxvgen and nitrogen for 
temperatures much above a few hundred degrees Kelvin, the Rankine— 
Hugoniot equations with constant specific heat ratio are not accurate for 
shock waves of Mach number 2 or greater. As the Mach number in- 
creases, this discrepancy also increases. For these stronger shock waves, 
the equilibrium static temperature behind the shock, 7’, , is so high that the 
vibrational modes of a diatomic gas become excited, thus changing the 
specific heat of the gas. (‘This static temperature of the gas is a measure 
of the random molecular velocity, irrespective of the mass flow as a whole.) 
If the shock strength becomes great enough, and 7’, high enough, dis- 
sociation of the molecules can take place, and, for still higher temperatures, 
even ionization may occur. In the present investigation, however, only 
vibration ettects are studied since 7, was low enough for dissociation and 
ionization ettects to be negligible. A method of calculating the shock wave 
properties, taking into account this variation of specific heat with tempera- 
ture, has been developed by Bethe & ‘Teller (1951). 

‘The process by which a sudden increase in energy is distributed among 
the various degrees of freedom of a gas is very complex. Such a sudden 
increase in energy occurs when a shock wave passes through a gas. How- 
ever, if the gas Is diatomic, the redistribution of energy is not instantaneous. 
The translational degrees of treedom arrive at a Maxwellian distribution in 
the space of a few collisions. It can also be shown that the energy of 
rotation approaches equilibrium with energy of translation very rapidly, 
again within several collisions. ‘This state is not the final equilibrium 
state of the gas, since, for molecules having vibrational states, the time 
necessary for the energy of vibration to come into equilibrium at the final 
temperature is known, in some cases, to be many thousand times that 
required for translation and rotation. ‘This time, which is characteristic 
of the process, may be represented by a relaxation time for vibration, 7. 

Pierce (1925), using the magnetostriction oscillator for determining 
sound velocities, was the first to measure sound dispersion. Herzfeld & 
Rice (1928) showed that dispersion could be explained by ascribing a time 
for readjustment of the vibrational energy which was long compared with 
the time for adjustment of translation and rotation. ‘This led to the 
designation of translation and rotation as ‘external’ degrees of freedom, 
and modes such as vibration or ionization, which take a relatively long time 
to adjust themselves to external conditions, as ‘ internal’ degrees of freedom. 
In the theory of sound dispersion the total specific heat C,, may be thought 


of as being made up of two parts: 


( = ¢ thi ( vib’ ( l ) 
where (’,’ is the specitic heat due to the external degrees of freedom, and 
C\, is the specific heat in vibration. Since the energy in the vibrations 


responds slowly to external changes, C\,,, depends on the frequency of the 
sound. ‘Thus, tor very low frequencies, C,,, 


can follow the etfects of the 
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compressions, but, for much higher frequencies, the time between com- 
pressions and rarefractions can become short compared with the adjust- 
ment time for vibration, and the effective specific heat of the gas does not 
include C,;,. This means that the effective value of y for the gas increases, 
which in turn means that the velocity of sound also _ increases. 
At intermediate frequencies, the energy in vibration partially follows the 
changes in pressure but is out of phase with it. This causes the sound 
absorption to pass through a maximum. ‘Thus, measurements of sound 
absorption as a function of frequency enable vibrational relaxation times to 
be determined. ‘The method is well described by Herzfeld (1953). 

A second method for measuring relaxation times was suggested by 
Kantrowitz (1942). In this method one allows a jet of gas to expand 
through a nozzle in a time which is long compared with the relaxation time, 
so that the full C, of the gas is active. If this jet is then stopped at the head 
of a pitot tube in a time which is short compared with the relaxation time, 
only the specific heat of the external degrees of freedom isimportant. ‘Chen 
the difference between the tank pressure before the expansion and the stag- 
nation pressure measured by the pitot tube is a measure of the energy lag in 
the gas. When the jet is stopped in a time of the order of 7, the relaxation 
time for this lag can be determined from the dimensions of the tube and the 
streamlines of the flow. Kantrowitz & Huber (1947) and Griffith (1950) 
have used this technique to determine the relaxation times for a number of 
gases, and the method seems well established. 

A third method for determining relaxation times, which makes use of a 
shock tube and interferometer, was suggested in 1950 by Griffith. 
Measurements of 7 for CO,, N,O, and methane have been made by this 
method at Princeton by Griffith, Brickl & Blackman (1956), and measure- 
ments on CO, and Cl, have been made by Smiley, Winkler & Slawsky (1954). 
The big advantage of this method over the two previous methods is that, by 
using various shock strengths, the value of 7 for widely different tempera- 
tures can be determined. ‘This serves to give a convenient check on a 
theory by Landau & Teller (1946) which predicts the dependence of 7 on 
the temperature. 


‘THEORY 

The Rankine—Hugoniot equations are exact expressions only for perfect 
gases with constant specific heat. Since the specific heat changes with 
temperature, it is clear that these equations are not accurate for strong 
shock waves. For argon, Resler, Lin & Kantrowitz (1950) have shown 
that this change in specific heat is slight up to .W, = 9; however, above VW, =9 
in argon, where temperatures greater than 8000° K are produced, ionization 
begins to occur and the specific heat may no longer be treated as a constant. 

For diatomic gases such as oxygen or nitrogen the dependence of C,, on 
temperature is of importance at about 500° K, and for carbon dioxide this 
variation of C,, is of considerable significance even at room temperature. 
‘The translational and rotational modes of O, and N, are fully excited to the 
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classical value of R 2 per degree of freedom per mole at room temperature, 
while vibration and especially dissociation and ionization are excited to 
such a small extent that the energy in these modes may be neglected. As 
the temperature of the gas increases, these internal degrees of freedom can 
become excited, with the vibrational mode being the first excited since it 
has the lowest energy. For diatomic gases the energy in vibration, £,,,, 
is given by 


exp(4T)—1? 


Evin (2) 
where = Av k is the characteristic temperature for vibration. For nitrogen 
() = 3336 K, while for oxygen @=2228 K. For oxygen C, increases, due to 
vibrations, by 1:27R for an increase of temperature from 300° to 3000° Kk, 
while for the same increase of temperature in nitrogen C’, increases by 
0-95R. In oxygen about 6”, of the molecules are dissociated at 3500°K, 
and in nitrogen only about 1”, are dissociated at 5000 K (according to Bethe 
& Teller 1951). Here, only the energy in vibration will be considered, 
since the change in C’, due to dissociation, even though appreciable at the 
strongest shocks obtained, is characterized by such a long relaxation time 
that these effects did not show up in the region behind the shock that could 
be observed. 

In the Bethe—Teller method for calculating the shock parameters for 
the case of variable specific heat, one considers a one-dimensional shock 
wave moving through a gas. ‘Then the conservation equations may be 
written : 


mass, =e} (3) 
momentum, Py + pity = Pot pots; (+) 
energy, Ey + py py 2= Ey +t} 2. (5) 


For most gases at ordinary densities the perfect gas law applies: 
p p=RT. (6) 


‘These four equations uniquely determine the equilibrium values of the 
shock wave parameters, independently of the manner in which the equili- 
brium is attained. 

From these equations, one is able to calculate the density, temperature, 
and pressure ratios as functions of .7, for any strength of shock, provided 
the dependence of (, on temperature is known. In order to facilitate 
these calculations, Bethe & Teller define a new quantity: 


Bip) =1+ Ex RT, 


where £,, is the total energy content per unit mass. If equilibrium exists, 
it can be shown that f is a function of temperature only. The energy 


equation can be rewritten in terms of 8: 


B, RT, +7? 2=B,RT, +73 2. (8) 
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‘The solution of these equations can now be obtained. The density ratio 
across the shock is given by 


Ty, (9) 
where 
b= B,— 1/2 —(B, — 1/2)7,/T,. (10) 
The Mach number is given by 
My By Ty) — (11) 


in which y, is the ratio of the specific heats at the temperature 7. For the 
calculations on oxygen and nitrogen, values tor 8 were taken from the 
reports on these gases by Woolley (1953a,b). It should be remembered 
that p, and 7’, refer to the final equilibrium values behind the shock wave. 
According to this viewpoint, the shock thickness includes the region that is 
necessary for vibrational equilibrium to be reached. ‘This transition 
region may be very broad and depends on the pressure and temperature 
behind the shock. 

Since the density and the speed of the shock are the two quantities that 
are measured directly by experiment, they are of the most importance as far 
as comparing theory with experiment is concerned. The pressure ratio 
across the shock may be obtained from the equation 


= pols (12) 
Figure 1 shows the pressure, density, and temperature ratios across a 
shock of speed M/, in oxygen, as calculated by the method outlined above. 
The same quantities, as found by the Rankine—Hugoniot relations with 
constant specific heat ratio, are plotted for comparison. ‘Table 1 gives 
numerical values tor the quantities used in the calculation, while table 2 
gives the same quantities for Ng. 

In figure 1, the curves for constant specific heat are labelled as p, p,, 
etc. State @ here represents the region immediately behind the shock 
which is characterized by equilibrium between translation and rotation 
but no change in the vibrational energy. For any gas, the variables at state 
a can be calculated from the Rankine—-Hugoniot equations if the contri- 
bution of the vibrational energy to C,, is omitted. For oxygen and nitro- 
gen, the Rankine~Hugoniot equations with y = 1-40 can be used to calculate 
variables in state a, since the amount of vibrational energy at room tempera- 
ture is negligible. 

The transition between state a and the final equilibrium is characterized 
by a changing energy of vibration, F,,,, for which Landau & Teller have 
derived the relaxation equation from general quantum considerations. 
‘hey assume that the molecules are free to vibrate in any quantum level, 
but that changes in energy occur for m= +1 only. If ky and ky, are the 
probabilities of dé-excitation and excitation per unit time respectively, 
hetween the ground and first excited states, then 

hoy = exp( —hv/RT) (13) 
F.M. E 
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from the principle of detailed balancing. Likewise, for other quantum 
states, 


= (14) 


and 


= exp(—Av RT). (15) 


9 10 


Figure 1. Density, pressure, and temperature ratios for shock waves in oxygen. 
The ratios with subscript a were calculated assuming constant Cp, while 
those with a subscript 2 are calculated by the Bethe—Teller method for variable 
Cp. 
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296 B, =3-500 


b Po! py M, Pal Pr 
300 3-500 | 0-040 | 1-034 | 1-019 | 1-033 | 
400 | 3-519 | 0-799 | 1-973 | | 1-956 | 402, 
600 3-593 613 3-373 2-434 256 ; 613 

SOO 3-680 2 076 4-230 3 134 3 075, 851 

1000 2-387 4-834 3-720 4-408 1073 

1200 3°853 2 613 5 273 4-698 1316 


1400 3-921 4-899 1555 


1600 3-981 2-926 5-883 5-148 5-049 1805 


1800 4-03 3-041 6-108 5-563 5-159 | 2060 


00 

to 
+ 


2000 4-082 3-138 6-300 5-948 5-253 2315 


2200 4-126 3°223 6:457 6-314 2555 


2400 +-168 3-298 6-615 6-063 2810 


4-208 3-366 6°749 6-999 5-444 3098 


2600 


Table 1. Shock parameters for oxygen. Values of B were taken from the report 
by Woolley (1953 b). 
Then, if y,, is the number of molecules in nth quantum state, 

dy, 
dt 
By multiplying equation (16) by mv, and summing over 7, one can 

obtain the expression for the rate of change of F,,,, with time: 
OF vib 


ot 


k,, in Yn oF k, +1,n +1¥n-° (16) 


= ky fl —exp(—hv/k (T) Evin} (17) 


where E£,,,, (7) is the energy in vibration when the gas is in equilibrium with 
the external temperature 7. If the external temperature does not change 
much during the process, the solution of equation (17) is given by 


Evin — Eyin(T) = C exp(—t/r) (18) 


with 


: =k,{1—exp(—hv/kT)}. (19) 
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For shock speeds otf W,=7:5 in oxygen and W,=10 in_ nitrogen, 
the assumption that 7’ remains constant in the above integration may be in 
error by as much as 20°,. 


T,=296°K = B, ~ 3-497 | 
| b | Po Pi | M, | Pa Py | 

| | 
| | | 
300 3-497 | 0-040 | 1-041 | 1-024 1-041 | 300 
| | | 
| | 
| 
400 | 3-500 | 0-782 1-938 1-547 1-938 | 401 
| | | | 
| > ? | 
| 600 3°521 | 1°542 3-237 | 2-397 = 604 
| 
| | | 
| 1000 3-624 2-237 | 4-539 3-624 4+:346 | 1074 
| | | 
| 
1400 | 3-746 | 2-614 | 5-267 | 4-592 | 4-850 | 1495 
| 
1800 | 3-852 | 2-859 | 5-747 | 5-418 | 5121 | 1980 
| | | | | | 
2200 | 3-938 | 3-035 | 6-092 | 6-151 | 5-295 2460 
| | | | | 
2000 | 4-008 | 3-167 | 6-351 | 6-815 | 5-415 | 2960 
| | | 
3000 4-065 3-268 | 6°551 | 7-425 | 5-500 | 3450 
| | | 
3500 | 4:123 | 3-365 | 6-742 | 8-130 | 5-575 | 4080 
| 
| 
4000 | 4:171 | 3-449 | 6-901 8-783 | 5-632 | 4690 
} | 
| | 
4500 | 4-211 | 3-513 | 7-035 9-392 | 5-675 | 5330 
| 5000 4:245 | 3-568 | 7-144 | 9-969 | 5-710 | 6010 


Table 2. Shock parameters for nitrogen. Values of B were taken from the report 
by Woolley (1953 a). 


Since the measurements that are carried out are actually on the density 
tield behind the shock, it is necessary to find an expression connecting the 
change of density with the relaxation time. 

If one assumes the transition region to be at a constant pressure, then 
the following relation between state @ and the final state holds: 

dE in = Cy, = — C,'aT. (20) 
For linear diatomic molecules, C,’=7R 2, and this may be considered as 
constant in the transition region, since translation and rotation are already 
fully excited. Equation (20) may then be integrated to give 
Ev» — Con —C, T+C,'T2, 


so that 


¢ 
| 
\ 
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Equations (17), (20), and (21) can then be combined and integrated to give 
T-T, 
T_T. =exp(—C,t/C, (22) 
From the perfect gas law, and the assumption that the pressure remains 
constant this may be converted into an expression for the density change : 


Ds 2 4 
Pe exp( = Pa exp(— Cc; | . (23) 


P2—Pa 


The term exp(—C,¢-C,’r) in equation (23) is about 0-1, or less, for 


a 


most of the experimental work reported here, so that to a first approximation, 
the whole term in square brackets may be set equal to 1. ‘Therefore, not 
only the temperature, but also the density, approaches the final equilibrium 
value in an exponential manner. 

The relaxation region is characterized by 7\,,, increasing from 7', to 
7,, and 7, decreasing to 7,. ‘The energy in vibration F,,,, also increases 
rapidly with time. ‘lhe density and pressure both increase to the final 
equilibrium value with the same time constant as temperature. ‘The 
change in pressure from state a to the final.value, as shown by figure 1, is 
quite small and may be neglected in many calculations. However, the 
relatively large change in temperature and density is evident in the figures. 
It is this large change in density that makes the interferometer particularly 
suited to these measurements. 

Landau & eller have also developed a theory that helps to give a 
better understanding of the relatively slow exchange by collisions of vibra- 
tional energy in comparison with rotational energy. According to their 
work, the effectiveness of a given collision in exciting or de-exciting a 
periodic degree of freedom depends on the degree to which the collision 
can be considered as adiabatic. In this sense, a collision is assumed to be 
strictly adiabatic if the degree of freedom under consideration undergoes 
an infinite number of cycles during the time of contact. In such a case, 
there is zero probability that the mode will change its state during the 
collision. However, if the collision is not adiabatic, that is, if there is only 
a finite number of oscillations during the perturbing effect of the collision, 
then energy may be transferred from the translation to the degree of free- 
dom being considered. ‘To express this fact mathematically, they use the 
ratio x= (24) 
where 7". is the duration of the collision, and 7;, the natural period of the 
degree of freedom. If x is of the order of unity or less, then only a few 
collisions are necessary to bring the degree of freedom into equilibrium, 
while if y is much larger than unity, many collisions will be needed for 
equilibrium. In their report, Bethe & ‘Teller show that x <1 for rotation, 
and therefore the rotational energy of a gas should come into equilibrium 
with the translatonal energy in just a few collisions. Experiments on 
the thickness of shocks in gases have shown this to be the case. 
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For a gas that has vibrational levels, this equilibrium between trans- 
lation and rotation is not the final equilibrium state, since the vibrational 
energy will also change. When one considers vibration, 7. in equation 
(24) may be written 


(25) 


where r is the range of strong molecular interaction, and 7’ is the relative 
velocity of the two colliding molecules. In equation (24), 7), is simply 
the reciprocal of 27 times the frequency of vibration, so that 


X= (26) 


For nitrogen at room temperatures this ratio is about 15. Since this is 
much larger than unity, it will take many hundreds of thousands of collisions 
to excite the vibrational states of nitrogen. 

When y >1, as is the case for vibration, Landau & ‘Teller show that ?,, 
the probability per collision of de-exciting a molecule in the first excited 
state, is given by 
Py. = Cexp(— x), (27) 


where ( is a factor which depends on the geometry ot the collision and ts 


arbitrarily set equal to 1 10, 


It may be noted from equations (26) and (27) that, as the temperature 
increases (7 increases), the probability of de-excitation per collision in- 


creases, which in turn shortens the relaxation time. 

Now equation (27) must be generalized to include all of the molecules 
in the gas by averaging over the Maxweli velocity distribution. When 
this is done, it can be shown that one arrives at the following expression : 

expf K, 723), (28) 
K, and A, are constants which depend on collision cross-sections and which 
cannot be accurately determined. Now P,, is related to k,, by the equation 
VP. =Ry, (29) 
where N is the number of collisions per unit time. ‘Therefore, 
| > 
- = .VP,,{1—exp(—hv RT)}, (30)) 
or since 
N pF 
T'®exp(K, —exp(—hAv kT}. (31) 

For relatively low temperatures (~< 1000°K) equation (31) shows that 
logr oc 7°18, and this affords a convenient method of checking the theory. 
However, as the temperature rises. the factor {1 —exp(—Av kT)! can 
become important and must be taken into account in trying to compare 


experiment with theory if the accuracy of the experiments warrants it. 
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In recent years Schwartz & Herzfeld (1954) have elaborated on this 
theory, and have arrived at an equation which allows numerical compu- 
tation of P,,. They give the result 


1 3 (hv \2 (RT\18 hy 


where « and «’ are related to molecular constants. 


EXPERIMENTAL EQUIPMENT 

In the course of normal operation, with helium at three atmospheres 
pressure in the high pressure chamber, and air at about 8 mm Hg pressure 
in the channel, experiments in the shock tube have been limited to ., =4 
in air. Experiments on vibraticnal relaxation times in CO,, methane, and 
especially NO have been verv successful with the present arrangement 
for two reasons. First, at the temperatures obtained behind shocks of 
this strength, there is sufficient vibrational excitation in these gases to 
cause a significant change in density from region a to the final equilibrium 
state in region 2. Second, the relaxation time is of such a value that the 
actual physical extension of the relaxing region behind the shock can be 
easily photographed and measured. 

Measurements with the Mach-Zehnder interferometer on shock waves 
in nitrogen showed that the density ratio across the shock corresponded to 
p; mstead of p, Hence, the vibrational adjustment has not taken 
place within the region of tlow behind the shock that can be studied by this 
means. ‘his was to be expected since it has been known for some time 
that relaxation times for diatomic gases are considerably longer than for 
polyatomic ones. An estimate of the lower limit of > for nitrogen showed 
that + -200usec at 900 K for atmospheric pressure behind the shock. 
Since the relaxation time decreases as the temperature or pressure increases 
(according to equation (31)), it was necessary to devise a method of obtaining 
stronger shock waves. It was not practical to do this by simply reducing 
the pressure at the downstream end below 8mm Hg, since the interfero- 
meter would have become increasingly insensitive. . 

After considering several possible methods for obtaining stroriger 
shocks, it was decided to construct a converging channel which could be 
mounted inside the present tube. ‘T'his plan permitted use of the inter- 
ferometer, together with the associated optical equipment and many of the 
existing electronic triggering circuits and devices. The incident shock, of 
speeds up to VW, =4, was formed in the channel of the tube and allowed to 
converge through an area ratio of 34:1 to a channel of constant area. 

The interferograms of the shock were taken after the wave had travelled 
about 20 inches along this final straight channel. Thus, if one assumes 
that the shock strength is constant fn the straight section, then, for a density 
ratio of 6 across the shock, there should be over 3 inches of uniform flow 
behind the shock by the time it reaches the centre ot the window. In 
practice, it was noticed that, as p, increased, irregularities in the How 
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became more evident in the photographs. ‘This is due to the increase in 
sensitivity of the interferometer as the density increases. However, in 
general, there were at least two inches of uniform flow behind the shock, 
and this distance was more than sufficient to measure relaxation times. 
‘The increase in shock strength was very considerable, and, for a wide range 
of initial pressure ratios across the diaphragm, the Mach number of the 
shock was almost doubled upon convergence. With helium as the driver 
gas, the upper limit for work in oxygen, nitrogen, and air was M,=7°5. 
For some of the work, hydrogen was used in the high pressure end in order 
to obtain shock strengths up to M,=10 in nitrogen. At these shock 
strengths, relaxation effects were evident in both oxygen and nitrogen. 

In order to check the Bethe—Teller method of calculation, two inde- 
pendent shock parameters had to be determined for each shock wave. The 
first of these measurements was the velocity of the shock, since from this 
velocity, and the speed of sound ahead of the shock, the Mach number can 
be determined. In order to measure the shock velocity, the time of passage 
of the shock past two points had to be determined. ‘The detecting device 
that was used consists of a very narrow strip of evaporated gold whose 
resistance changed when it was heated up by the passage of a shock.* 
‘These gold films were of the order of 100 to 300 atoms thick, 1 2mm wide, 
and 7mm long. ‘The resistarice of the films was between 15 and 20 ohms, 
although films of widely varying resistance could be made to work. In 
use, a 22 volt battery and 2000 ohm resistance were set up in series with a 
tilm, then, when a shock passed over the film, a small change in voltage 
occurred across it. ‘The rise time of this voltage pulse is determined only 
by the transit time of the shock across a film. As a test of this, a shock was 
allowed to travel lengthwise along a film one inch long, and the signal was 
put on an oscilloscope (see plate 2). “The velocity of the shock, determined 
from the transit time of the shock along the film and the length of the film, 
agreed with the velocity measured by other means to within 1°,. The 
narrow gold tilms that were used in the channel were oriented perpendicular 
to the motion of the shock, and the signal rise time was about 0-5 ;:sec or less. 

‘The gold tilms were made by evaporating a thin laver of gold directly on 
to sound recorder mending tape. ‘This tape was chosen because it is the 
thinnest of the commercial adhesive tapes, with a thickness of about 2 mils. 
Strips of this tape could then be cut to the proper width and pressed in place 
gently, without harming the gold surface. The thickness of the films was 
estimated by weighing the filament both before and after the gold was 
evaporated and measuring the distance from the surface to the gold. From 
this information and the assumption that every atom of gold that strikes 
the surface sticks to it, a simple calculation gives the thickness of the film. 
It was found that the gold strips were extremely durable, showing no_ 
measurable change in resistance from shot to shot, and would last for dozens 
of shots if they were not scratched by a piece of the diaphragm. These 
strips of tape with the evaporated gold on the surface were stuck fast to 


* The development of this technique was aided by information trom R. J. Emrich, 
of Lehigh University, concerning a heat flux meter made of thin gold foil. 
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lucite plugs mounted flush with the surface of the upper wedge, as shown 
in plate 1. Contact from the film to the lead through screws was made by 
painting DuPont silver paint over the ends of the strip and the screw. 
This arrangement proved very satisfactory and enabled a broken film to be 
replaced with a minimum loss of time. 

Provisions were made so that five of these thin films could be mounted 
along the upper wedge, although in most cases only two films were used. 
‘These were mounted five inches apart and situated so that 7, was measured 
just before the photograph of the shock was taken. 

The voltage output from the thin films amounted to 10 mv or more for 
the shock strengths used. ‘This was more than 100 times the noise level of 
the apparatus. A typical oscilloscope trace of a thin film pulse is shown in 
plate 2. Since the signals were flat topped, they were differentiated before 
entering the chronograph. ‘The timing traces on which the signals appeared 
had 5 usec markers and the time of transit of the shock between the 2 films i 
could be read to 0-2 sec, giving an error in the time measurement of about 
+0-4°, for M,=7. ‘The distance between the films was measured to 
within +0-2°,, so that the error in the Mach number determined in this 
manner is about +0-6%, at .W@,=7. As the Mach number decreases this 
error decreases, since the time measurement becomes increasingly accurate. 3 
For many shocks two separate measurements of the shock velocity were 
made, and these two measurements usually agreed very well. 

‘The second measurement that was taken on each shock was tie measure- 
ment of the density of the flow, which was made photographically with the 
use of the Mach-Zehnder interferometer. ‘lhe illumination for the inter- 
ferograms was obtained from a barium titanate spark with an effective 


duration of less than 0-3 ;.sec, which was sufficiently short to stop the motion ¥ 
of the shock. When charged to 10,000 volts, the illumination was intense . 


enough for all purposes. ‘The electrodes were made from stainless steel 
rods which were cleaned from time to time in order to keep the sparking 
reliable. ‘The time at which the spark fired could be varied by an adjust- 
able delay, the delay being triggered by the signal from one of the thin films. 

From the interferograms the density in the region behind the shock 
can be determined from the fringe shift in the given region, as described 
by Bleakner, Weimer & Fletcher (1949). If 5 is the shift of a fringe from 
its position at a given point in the field, and p, is the pressure ahead of the 
shock, then the density p at this point is given by 


= —8 (33) 

Pi Pi 
where W is a constant depending on the gas. The value of W can be 
either determined experimentally, by counting the fringe shift 6. past a 
given point as the pressure in the shock tube is slowly changed and using 
Py , or calculated from the equation 

Un, 1) 


the expression 


Ww (34) 
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where / is the width of the tube, ,,7,, and p, are the index of refraction, 
temperature, and pressure of the gas at S.T.P., and A is the wavelength of 
light used. Values of HW’ obtained by these two methods differ by less than 
EXPERIMENTAL RESULTS 

Shock waves in argon, oxygen, nitrogen, and oxygen nitrogen mixtures 
have been studied by the methods just described. ‘Typical interferograms 
of shock waves in these gases, and the flow behind the shocks are shown in 
plates 3, 4, and 5. From these pictures and the simultaneous measure- 
ment of z7,, a check on the Bethe-'leller method of calculation has been 
carried out, and the vibrational relaxation times have been measured. ‘The 
density ratios p, p, and py p, were determined from the pictures and the use 
of equation (33) where 6, is the fringe shift from in front of the shock to 
region « and 0, is the tringe shift to the equilibrium region. — The fractional 
fringe shift can be measured, reproduceable to + 0-04 fringes, by the use of 
an optical comparator. However, the determination of 0, and 0, are not 
always this accurate. ‘The value of 0, depends on the point at which the 


fringe meets the shock. “This measurement is greatly atfected by the 
8 Prism 
Z , Window 
—s 
|) 
_|_LOPTICAL AXIS 
Shock Front | | 
} 
Shock TubeWall”! | 8 Prism 
Figure 2. A schematic diagram of the offset prism technique which shows how the 
light is sent diagonally through the shock. “Phe interferometer has been 


ymitted., 


apparent width of the shock in the interferogram (due to the spark duration 
and any slight angle the shock may have with the optical axis) and also by 
the rapid change of fringe position immediately behind the shock for short 
relaxation times. ‘hus the value of 6, mav be in error by as much as 0:2 
fringes. ‘The value of 5, for a gas showing relaxation may also be in error 
due to curvature of the undisturbed interferometer fringes. Since 6, must 
be measured | inch or more behind the shock in some cases, errors in 6, may 
be as high as 0-1 fringes. In general, errors in 6, increase as the relaxation 
time decreases, and errors in 6, increase as the relaxation time increases. 
In order to identity the integral number of the fringe shift across the 
shock, a technique using small angle offset prisms was used to send light 
obliquely through the shock.* A schematic diagram of the offset prism 
arrangement is shown in figure 2. This method of identifying the total 
fringe shift is a slight modification of a technique using an off-axis 


“The offset prisms were first used tn this laboratory by D. R. White. 
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light source. Plate 4 shows an interferogram in which the shock can be 
seen in the channel and, at the same time, in the prism. It should be noted 
that a fringe can now be traced through the shock and the total fringe shift 
established. 

As a check on the validity of the method of measuring M, and p/p), 
measurements were first made on a number of shocks in argon up to M@,=8. 
Since there is no excitation of any internal degrees of freedom (y= constant) 
for these shock strengths, the values of p,p, as a function of M, can be 
calculated from the Rankine-Hugoniot equations with essentially no error. 
‘The measured values of p, p, and VM, were compared with the calculated 
values, and the agreement proved satisfactory. A shock wave passing 
through argon is shown in plate 3. It can be seen that the fringes do not 
shift in the region behind the shock, indicating that the density is constant. 
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Figure 3. The experimental and theoretical density ratios in oxygen. 


A comparison of the two interferograms of shocks in oxygen with the 
one in argon (plate 3) shows a considerable difference in the region behind 
the shock. In oxygen, the fringes shift from the value at the shock to a 
higher final value in the equilibrium region. This shift in fringes reflects 
the density change in the transition region from state a to the final equili- 
brium state. From the interferograms the values of p, p, and py p, have 
been determined for a number of shock waves in oxygen. - 

Values of Ap p, where, Ap is the difference between the experimental 
and theoretical values of p, for the one case, and the difference between the 
experimental value of p, and the theoretical value of p,, for the other case, 
are plotted against the experimentally determined Mach number in figure 3. 
The difference between the theoretical values of ps/p, and p,/p, has been 
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plotted for comparison. It can be seen from the figure that the experi- 
mental values of p, p, agree very well with the theoretical calculation. 
However, the results also show that the measured values of p,/p, lie con- 
sistently below the calculated values by about 1 to 2°. Since there is only 
one vibrational mode in oxygen, and since any dissociation effects would 
tend to make the density higher instead of lower, this discrepancy in pop, is 
believed to be due to impurities and errors in measurement. The purity 
of the oxygen was not very high, and, for the strongest shocks, there was 
about 5°, of nitrogen present in the gas ahead of the shock. As the value 
of p, increased, this percentage of nitrogen decreased. It is difficult to 
obtain a tube of oxygen of high purity with the present arrangement, since 
the pumping system takes the pressure down to only about 0-5 mm Hg, and 
the oxygen cannot be flushed through the shock tube like argon or nitrogen, 
since it would react with the pump oil. Hence, the purity of the oxygen 
varied between 95", and 99°, for the work reported here, with the main 
impurity being nitrogen. ‘The nitrogen still has essentially no energy in 
vibration by the time the density ratio reaches p,/p,, so that the role of nitro- 
gen as an impurity should be to cause the density in the region after the 
oxygen relaxes to be lower than the calculated value. Therefore, this may 
explain part of the discrepancy, but it is believed that another error was 
present in the determination of the pressure which may have been as large 
as 1”, for some cases. 

The relaxation time in oxygen has also been determined from the inter- 
ferograms. Equations (23) and (33) can be combined to give an expression 
relating the fringe shift in the transition region to the relaxationtime. ‘Thus, 


oO oO ) 


-exp( 1+ — exp( — C,,t/,C7) (35) 


Os 


Since the denominator is nearly equal to 1, by taking the natural logarithm 
one obtains the result 


log d= —C,t/Cir+K (36) 


or, since t= .x 7, where x is the distance behind the shock and @ 1s the average 
velocity of flow behind the shock, 


vlog(o,/6,) 


(Xp, (3 7) 


This equation may now be used to calculate + from the interferograms. 
Figure + shows some experimental plots of logd versus the distance 
behind the shock front. ‘The fact that these points lie approximately on a 
straight line justifies the replacement of the term in square brackets in (35) 
by unity. These plots give one confidence in this method of measuring 7. 
The values of (x, —.x,) and log(6,/5,) are determined from plots such as 
these, and @ is taken as the arithmetical average of 7, and v9. 
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Figure 4. Typical plots of fringe shift versus distance behind the shock. § and x 
are in arbitrary units. 


The values of + for oxygen using this method have been calculated for 
a number of shock strengths and the logarithm of the results has been 
plotted against 7-13 in figure 5. To the extent that the factor 
{1—exp(—Av/kT)} in equation (31) can be neglected, this sort of plot gives 
a direct check on the validity of the Landau—Teller theory. The relaxation 
time calculated from equation (37) is 7 at the pressure py, so that these 
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values must then be corrected to atmospheric pressure for the plot in 
figure 1. The temperature 7, which is used for this graph is the average 
of 7, and 7,. It can be seen from the figure that the experimental points 
appear to lie on a slightly curved line instead of a straight line. 
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Figure 5, Vibrational relaxation times for atmospheric pressure in oxygen. 


‘The oxygen, and also the other gases, were dried in a 4-foot long drying 
system containing CaSO,. Since the percentage of water vapour de- 
creases as the pressure of the gas in the drier increases, the gas was flowed 
through this pipe at 20 to 30p.s.i. However, with the large quantity of 
gas to be handled in the shock tube, it is difficult to obtain the same purity 
as one can in a small volume. For the experiments in oxygen, the amount 
of water vapour varied trom about 1 part in 4000 for WW, =7, p, = 10mm Hg, 
to about 1 part in 20,000 for M7, =4:25, p,=50mm Hg. Since it is known 
that water vapour is very effective in bringing about vibrational equilibrium 
in many gases, it was believed that, for the stronger shocks, the greater 
amount of water vapour was lowering the relaxation time. Hewever, data 
on the relaxation times in the air, in which the gas was not dried, seem to 
show that, for strong shocks, water vapour in these concentrations does not 
shorten z+ by a measurable amount. Later experiments also showed that 
nitrogen did not have a large effect in reducing the relaxation time of the 
oxygen, so that the effect of impur.ties on 7 is believed to be very small. 
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Results of measurements of p,/p; in the range of M,=5-10 for N, are 
plotted in figure 6. For the pictures in which relaxation was evident, 
po/p; Was measured, and these values are also plotted in figure 6. As in 
figure 3, the values actually plotted are \p/p,, where Ap is the difference in 
p, as explained before. The agreement between theory and experiment 
is within the experimental limits for these values. It can be seen that, for 


the stronger shocks and smaller p,, the accuracy of the ex »« » ents decreases. 
15 
Experime: entc! - in. 
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Figure 6. ‘The experimental and theoretical density ratios in nitrogen. 


The purity of the nitrogen used in these experiments was considerably 
greater than the purity of the oxygen. ‘The shock tube was pumped out to 
0-5 mm Hg pressure and then nitrogen which had been dried by the CaSO, 
drier was flushed though the tube for about three minutes before the final 
value of p, was set. This technique should keep the presence of water 
vapour in the experiments to less than 1 part in 10°. 

Relaxation times in nitrogen have also been calculated from pictures 
such as in plate 4 by the method outlined above. The results of these 
measurements are shown in figure 7, in which log7 is again plotted against 
7-13, as an experimental check on the Landau—Teller theory. From the 
few points plotted, it appears that the data check the theory within the 
accuracy of the measurements. 

Figure 8 shows a much wider range of z and 7, and includes relaxation 
times from shock tube measurements in oxygen, nitrogen, and CO,, Pitot 
tube measurements in nitrogen, and sound dispersion measurements for 
oxygen (Kneser & Knudsen 1935, and Knotzel & Knotzel 1948). 
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Figure 7. Vibrational relaxatron times for atmospheric pressure in nitrogen. 


From the relaxation times, which are directly calculated from the 


experimental data, one can also obtain transition probabilities. In the 
theoretical work by Herzfeld & Schwartz (1954) on oxygen and nitrogen, 
theoretical values of P,, for temperatures up to 3000° K are reported. 
Consequently, for the purpose of comparison, it is convenient to express 
the results in terms of Pj. as well as 7. 

From equation (30), 


= NP,,{1—exp(—Av/kT)}. 
But 
N=N (38) 
where.f is the number of molecules per unit volume, q is the collision 
cross-section, and @’, the mean relative velocity of the molecules, is 


k (39) 


m 


Now values of collision cross-sections are not well known and depend 
greatly on the type of molecular model that is used. It is also known that 
the value changes with temperature, decreasing to approximately a constant 
value as the temperature increases. Therefore, it will be assumed here 
that, in the temperature range under consideration, the value of g is constant. 
Values of .V as a function of temperature were calculated assuming that, for 
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Vernon Blackman, Vibrational relaxation in oxygen and nitrogen, Plate 2. 


Plate 2. T'ypical voltage traces from the thin gold films. The sweep speed of the | 


scope is 10 ye sec cm. 
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Vernon Blackman, Vibrational relaxation in oxygen and nitrogen, Plate 3. 
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Vernon Blackman, Vibrational relaxation in oxygen and nitrogen, Plate 4. 
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Plate 4. Disturbed and undisturbed interferograms showing total fringe shift 
and relaxation in N, 
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Vernon Blackman, Vibrational relaxation in oxygen and nitrogen, Plate 5. 
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collisions between oxygen molecules, g= 3-6 x 10° cm?, and, for collisions 
between nitrogen molecules, g=4:1 10°%em*. From these values of N 
and the average experimental values of 7, P,, was calculated for a number of 
temperatures. These values are tabulated along with those of Herzfeld & 
Schwartz in table 3. Values of P,, from the Kantrowitz & Huber (1947) 
measurements on nitrogen, and those of Knétzel & Knotzel (1948) in oxygen, 
are also included in the table. ‘The agreement is really excellent considering 
previous attempts to calculate relaxation times. 
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| Pyo exptl. | Pp calc. 
Oxygen cm* 
| 288 
900 1-1x10-> 
1200 2°4x 10-5 
1800 9-8 «10-5 8-6» 10-5 
2400 x10 5-5 x10— 
3000 
Nitrogen g=4:1X<10- cm? 
600 3x10 17 x 10-!° 
1800 
3000 1-6 «10-5 
4000 9-7 x 10-5 
5000 25x10" 


Table 3. Experimental and theoretical values of Py. q is the gas-kinetic collision 
cross-section. ‘The experimental values for oxygen are due to Knétzel & 
Knétzel (1948) ; those for nitrogen are due to Kantrowitz & Huber (1947) 
The calculated values are taken from Herzfeld & Schwartz (1954). 


Investigations were also conducted on oxygen-nitrogen mixtures, and 
the effect of N,—O, collisions on the approach to equilibrium for O, was 
determined for a narrow range of temperature. ‘lhe oxygen and nitrogen 
for these experiments were mixed in a tank outside the shock tube and 
allowed to stand for several hours. By this means, it is believed that the 
gases were adequately mixed before being used in the tube. Values of 
ps p, for these mixtures for the intermediate state have been calculated by 
the Bethe—Teller method. It is assumed that, in this intermediate state, 
the energy in the vibrational states of nitrogen had not had time to change, 
so that the value of £.(N,) is still 3-50, which is the value appropriate to 
300°K. The value of 8 for the mixture is then given by 


p=> D:Pis (40) 


where p, and £, are the partial pressure and energy content of the 7th com- 
ponent. ‘The values of p, p, for this intermediate state have been measured 
for four different mixtures of oxygen and nitrogen including air. ‘These 
experimental values, together with the theoretical curves, are plotted in 
figure 9. Values of p, p, have also been measured for these mixtures and are 
included in the figure. ‘The numbers plotted here are the actual values of 
p p, instead of the changes in this quantity. 

The relaxation times of the oxygen in the mixtures have been calculated 
fora number of pictures like those in plate5. As a first approximation, these 
relaxation times were reduced to atmospheric relaxation times by assuming 
that O.—N, collisions were not effective in bringing about equilibrium 
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Thus, it was assumed that the relaxation was taking 


place at a pressure equal to the partial pressure of the oxygen. ‘The value 
of 7, corrected for this pressure, is too low, which means that the relaxation 


takes place at a higher effective pressure. 


Now, if one assumes that O, — N, 


collisions are 40", as efficient as O, — O, collisions in bringing about equili- 
brium in the oxygen, then the relaxation time, as shown by the experimental 
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Figure 9. The density ratios in O,—N. mixtures. 


points in figure 10, is in good agreement with the average curve for oxygen. 
For higher efficiencies for O,— N, collisions, the value of + becomes cor- 
Therefore, for this range of temperature at least, 
O, — N, collisions are only about 40%, as effective in transferring a quantum 


respondingly higher. 
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of energy to the oxygen as O,—O, collisions. Previous estimates have 
placed the effectiveness of O, — N, collisions at about 5 times that of O, — O, 
collisions. 

Most of the experiments were carried out on dried gases. However, 
some experiments were conducted in undried air to see if the effect of 
water vapour could be detected. ‘The relaxation times for 3 different 
shock strengths in undried air are also included in figure 10, and it may be 
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Figure 10. Relaxation times in O, —N, mixtures for atmospheric pressure assuming 
O.—N, collisions are only 40%) as effective in transferring energy as O,—O, 
collisions The curve is the average value of 7 for pure Oy. 


seen that the water vapour did not have a measurable effect on 7. The 
amount of water vapour in these shots was estimated to be 1 to 3 parts in a 
thousand. It is known that the effectiveness of impurities in bringing 
about equilibrium decreases as the temperature increases, and it is believed 
that this is the reason the nitrogen and water vapour had a smaller effect on 
the oxygen than anticipated. For example, if an O,—H,O collision is only 
50 to 100 times as effective as an O,—O, collision, this would not affect 7 
enough to be detected in these experiments. 
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SUMMARY 

In this paper the problem is considered of long gravity waves 
approaching a semi-infinite barrier which extends parallel to the 
wave crests, the whole system being in rotation. It is well known 
that, when the rotation is zero, there is a ‘shadow region’ behind 
the barrier in which the disturbance diminishes rapidly with 
distance from the edge. However, it is shown that the rotation 
gives rise to an additional wave in the shadow region. ‘The crests 
of this wave are at right-angles to the incident wave, and it travels 
along the barrier without attenuation in that direction. ‘The 
amplitude falls off exponentially with distance from the barrier, as 
in a Kelvin wave. The amplitude at the barrier may exceed 
that of the incident waves. 

The problem arises in connexion with the propagation of tides 
and storm surges in the ocean. 


1. INTRODUCTION 

The work described in this paper arises from certain aspects of an 
investigation into the origin of storm surges. A particular example of 
such a surge is the raising in sea-level which occurred progressively round 
the North Sea and produced such serious flooding in February 1953, but 
it is known that many surges pass unnoticed because they do not result in 
serious destruction (Corkan 1948, Rossiter 1954). Some may be gener- 
ated by atmospheric disturbances over the North Sea itself, but others 
appear to originate outside the sea and then propagate into and around it 
in an anti-clockwise direction as free waves on the rotating earth. 

Here we are concerned with a possible mechanism for the propagation 
of these free waves from a region to the west and north-west of the British 
Isles into the North Sea. However, the results may have more general 
applications in oceanography and meteorology. 

The particle velocity associated with long waves on a rotating earth is 
itself rotatory when conditions are uniform along the wave crests (Proudman 
1953, p. 262). If sucha system of simple harmonic plane waves is incident 
normally on a semi-infinite barrier, one might expect that after the waves 


have passed the barrier their transverse velocity components act as a 
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source for waves propagating into the region behind the barrier. In the 
rather similar problem in acoustics, there is no such source and this region 
is ‘shadowed’ from the incident waves, the only disturbance being due to 
diffraction effects arising at the edge of the barrier and dying out with 
distance away from it (Lamb 1932, p. 538). 

In this paper, we investigate the effect of the rotation on the disturbance 
behind the barrier. It is found that there is indeed a system of progressive 
waves travelling along the barrier when it lies in the right half of the plane, 
as in figure 1 (the sense of rotation is assumed to be anti-clockwise as in 
the northern hemisphere). The crest height of these secondary waves is 
not uniform but decreases exponentially away from the barrier as in waves 
of Kelvin type (Proudman 1953, p. 253). The amplitude of the secondary 
waves is not attenuated with distance along the barrier. We calculate the 
amplitude of the waves at infinity, and show that for a certain range of 
frequencies they may exceed the amplitude of the incident waves. The 
disturbance along a path extending trom the edge of the barrier in any 
other direction dies out rapidly with distance. 


2. ‘THE EQUATIONS OF MOTION 

In this paper, we assume, as in most long-wave theory, that the vertical 
acceleration is small (Lamb 1932, p. 254). The equation of motion in the 
vertical direction then reduces to the hydrostatic equation p=gp,(z +), 
where p is the pressure, z is the depth below the mean free surface, ¢ is the 
elevation of the free surface above its mean level, and o, is the density. 

The horizontal equations of motion become, on substitution of the 
hydrostatic equation and neglect of the non-linear and viscous terms 
(Proudman 1953, p. 220), 


ou 

VU= 

ot Ox | 

(1) 
ov | 

a” dy’ 


where u, 7 are the components of velocity in the horizontal (x, y)-plane, and 

are functions of x, y, t only. The Coriolis parameter y is equal to 2w sind, 

where w is the angular velocity of the earth and 4 is the north latitude. 
The equation of continuity is (Proudman 1953, p. 220) 


Ou av 1 
(2) 


where h is the mean depth of water, which is assumed constant and large 
compared with 
From these equations, a differential equation for ¢ can be derived: 
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If ¢ is assumed to have a simple harmonic time factor exp( —/ct), equation 


(3) becomes 


sa 2¢=(), (4) 
where k? = (0? — y?) gh, 
and the elevation ¢ is now to be taken as a function of x, v only. Through- 
out the paper we assume o> y. 
090 


incident wave 


Figure 1. 


The conditions of the problem are that waves of the form 
aexp[i(kv—ot)] are incident, from y= — on a semi-infinite barrier 
placed at v=0, x —-0 (see figure 1). The boundary condition on the 
barrier is 7=. ‘lo satisfy the radiation condition at infinity, it is easiest 
to use Copson’s condition that k has a small positive imaginary part, and 
that the secondary waves due to the presence of the barrier tend to zero as 
r-» «x where r=(x?+.1°?)!? (Baker & Copson 1950, p. 154). When the 
solution is completed, we may: let the imaginary part of k tend to zero. 
This, or an equivalent condition, is usually assumed in practice, although 
it appears never to have been rigorously justified in the case of a boundary 
extending to infinity (Peters & Stexer 1954). 

By substitution of the boundary condition on 7 into equations (1), the 
boundary condition on ¢ is found to be 


5 
(3) 


a OX 


where p=yja<l. 
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3. THE INTEGRAL EQUATION 
We follow the method used by Karp (1950) and Copson (1946) in 
deriving an integral equation which, when solved, enables ¢ to be found. 
The elevation ((x,\) in the interior of a contour. C is first expressed in 
terms of a Green’s function and the boundary values of ¢(«,y) and 
0¢(x,v),en, where 6 én denotes differentiation along the outward normal to 
the boundary. ‘Thus, 


5 


C(x, V) = | 5 Xoo) — G(xyy; Xoo) > dsy, (6) 
0 0 4 


where ds, is the element of arc length and G(x,v; xp, yp) satisfies 
eG 


Oy? 


+ —d(x—x9)8(y — Vo), (7) 


6(x) being the Dirac delta function. ‘lhe Green’s function thus represents 
the effect of a source at v= Xp, V= Vy. 
We choose the free-space Green’s function 
XorVo) (x = Xo)? (v = Ver (8) 
where H’ is the Hankel function of the first kind (Watson 1944, p. 73). 
The contour C is a large circle about the origin whose radius ultimately 
tends to infinity, but it is indented along the positive x-axis to exclude the 
barrier (figure 1). "Then 
C(x,y) = = -G (cl) dx, + — —C—} ds, (9) 
| u,- 0 } On, On, 


#0 


where [¢] = ¢..(v,0)—¢ (x,0), the difference between ¢ on the two sides of the 
barrier, and (’, is the circular part of (. If we now take ¢’ to represent 
the secondary waves arising from the incidence of the primary waves 
aexp(tky) on the barrier, equation (9) may be written as 


(x,y) = | [¢] dx,+| ds, + a exp(thy). 
J OV, Jar Ome On 

Finally, by the radiation condition, the integral round C, tends to zero as 

the radius of the circle tends to ~, and 


AG 
C(x,v) = | (Ic dx, +a exp(tky). (10) 
4 


~ 


The boundary condition on ¢ (and therefore on [¢]) is now substituted 
from equation (5) into equation (10); and, after am integration by parts, 1 
follows that 


0G 
| 0 ? Yo=V ) (1 !) 


if it is assumed that [¢] is bounded and that it is zero at the origin. This 


imposition of a specific behaviour on ¢ as x — 0 is necessary in this type of 
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problem, and a condition of continuity at the origin would seem to be the 
most appropriate in the present case (Karp 1950). 
When the operator (¢ dy—ipé dx) is applied to equation (11) on y=0, 
it follows from the boundary condition, equation (5), that 


which is valid only for x0. 


4+. "THE SOLUTION OF THE INTEGRAL EQUATION 
Equation (12) is solved by the Wiener—Hopf method (Titchmarsh 1937, 
p. 339), of which only the outline is given below. 


Let 
f(x) =0 
| 
g(x)=a function defined by (14) | 
—() x >0, f 
(13) 
q(x) =0 0, 
= ak x >0,} 
(fa a 
I(x-—x,)=<(— —ip— ]G> 
Then equation (12) may be rewritten as 
| f (x9) Ux — x9) , (14) 


which now holds for all x, and defines g(x) for x <0. 
Equation (14) may be solved by taking the Fourier transform of both 
sides of the equation. We then find that 
G(x) = + L(x) F(«), (15) 
where G(x), O{z), L(x), are the Fourier transforms of g(x), g(x), U(x), 
f(x) respectively. Now, 
ak 
Q(z)=| g(x)exp(—tux) dx = atkexp(—tux) dx =—, 
/0 a 
ind Q(«) is regular for .7}%}<0. Similarly, we have (Erdelyi 1955) 
k* — 47(1 — p?)} 


2(k? —«?)12 


L(a)= 


> 


and L(x) is regular for |.4{x}]<./{k}. Next, we assume that f(x) = O(e'*), 
with .4{k,' >0 as v ++ x, so that f(x)is bounded at + x. It then follows 
that F(x) is regular and bounded for .4{x}<0. Similarly, it follows from 
equation (14) that | g(v)|=Of{ exp(—-/ {k} |x| )} as x — 0, and therefore that 
G(«) is regular and bounded for 4 {a} > —.4{k}. The domains of regularity 
of the various transforms in the complex x-plane are shown in figure 2. 
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Equation (15) now becomes 


1{k? — — p?)} 
— F(a). (16) 


G(«)= + 


The object now is to factorize L(x) into two parts L(x) and L_(«), one 
regular in an upper half-plane of the complex variable «, the other in the 
lower half-plane, and then to express equation (16) as the equality of two 
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Figure 2. Domains of regularity of the transforms. 


functions that are regular in two half-planes which have a common strip of 
regularity. The resulting equation is 


(17) 


The left-hand side of equation (17) is regular for 7 |x} >—/}hk}, and the 
right-hand side for .4‘z}<0. Therefore, as both sides are regular in the 
strip 0S {x'> —J/‘k}, they define a function E(x) which is regular over 
the entire z-plane by analytic continuation. 

It follows from equation (17) that, as G(x), F(x) are bounded in their 
respective half planes, E(x) is O(|~|+*) for the lower half plane and O(|z[-¥?) 
for the upper half plane. Thus E(x) is O(|z|!*) as It then follows, 
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from an extension to Liouville’s theorem, that E(«) is a polynomial of order 
less than 1/2 “of therefore a constant. This constant is zero because 
= as |x| —- oc in the upper half plane (Karp 1950). Therefore, 


Wk—a) F(a) + 
— 
so that F(a) = (18) 


— 
Although the explicit form of f(x) will not be required, it could now be 
calculated by means of the Fourier reciprocal theorem : 


fix) = | F(a)exp(iax)dz O<e<JS jk}. (19) 


‘THE SOLUTION C(x, y) 


From equation (11), we have 


0G 
(x,y) = | + ip 
0 


) - +aexp(iky), 


Vo CX 
= | f(x») m(x— xy) dxy + aexp(tky), (20) 
where m(x) = ts + y?)!2], (21) 
and its transform is 
= + 2) exp! i| — }, (22) 


which is regular in the strip |7{x!}]-~.7{k!. The Fourier convolution 
theorem states that 


F(a) exp(txx) dx. 


By this inane equation (20) may be written as 
(x, Vv) = 5-| 
where s* = 1--p*. For convenience we now split the right-hand side of 
equation (23) into three parts: 


where 

=aexp(thy), 

2a ; 4)! 2(k- 45) 


+ | v|(k?— dx. 
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We consider first the expression ¢,. Copson (1946) has shown how 
integrals of this nature may be transformed into more recognizable forms. 
Let v=pcos@, y=psin@, and consider the integral in the first quadrant 
0 <@-7/2. Then the line of integration may be deformed to that given by 
=-0 
/ 
/ 
/ 
T=+0 
= 00-18. 
00 


Figure 3. Paths of integration in a-plane. 


x=kcos(6+ir), with z real and with the limits — ~, x (figure 3). After 

some manipulation it is found that 

| ah sin }(@—78)cosh} 
0 

sin }(@—7B)cosh $7 

» cosh 7—cos(@—7f) 

exp[7kp|cosh 7 — cos(@—78)}] dr, (26) 


xp(ikp cosh 7) dr, (25) 


_ 2ia 2ia 


sinh $8 exp|ikp cos(@—1B)} | 


where cosh B= 1)s. 
Consider now the integral 


Then differentiation under the integral sign leads to 


= “icosh }rsin exp[ix(cosh 7—cos 4)] dr, 


since this integral is uniformly convergent for v>6, where 6 is any small 
positive number. Hence 


dI(x , exp(27x sin? } 
= = (5) texp(tiz) sin P( id) 


~cosh3rsin 


9 cosh cos dr. (27) 
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This may be integrated to yield 


-(2a)# sinh $id 
I(x) =7'2exp(jiz) | exp( —7A®) dA + constant. (28) 
~O 


On putting x=0 in equation (27), we find 
1(0) = 

Thus, as /(x) given in equation (28) is uniformly convergent for » > 0), 
(2x) sinh hid 


I(x) =7'* exp({iz) | exp(—1A?) dA + 


-(22)* sinh sid 
= exp(jiz)| exp( — 7A?) da. (29) 


- x 


On substituting equation (29) into equation (26), we obtain 


— ¢\ 1/2 
exp[— + tkp cos(@ — 7)] 


-(2kp)? sinh 4(8 10) 


dA, (30) 
or 


“2 
=(— ‘) exp(— fiz + explis" | exp(— A?) an], (31) 


where z= (2kp)!” sinh +76). 

When considering ¢, in the other quadrants, the same type of trans- 
formation is used except that, for x <0, the z-contour is doubled back round 
the branch line which radiates from x= —k. ‘The results may be ex- 
pressed by the same equation (31) provided @ is allowed to range from 0 to 
27. 

The transformation of ¢, proceeds along similar lines. Neither ¢, nor 
¢, contains the rotational parameters at ali. ‘These functions may be 
expected to represent the usual diffraction and reflection effects of acoustics 
when the boundary condition is that the normal gradient of the dependent 
variable is zero on the barrier. 

Before stating the solution, it should be mentioned that the deformation 
of the contour when x~-0 involves crossing a pole at <=0. This pole 
contributes a term which cancels with ¢, for y>-0 (in the shadow zone) and 
reinforces it for y <0 (total reflection). When the pole is not crossed. 

The solution ¢, + ¢, may now be stated in the form 


a | 
+ — exp(— liz), exp(tky) 


-(2kp)$ sin + 0 -(2kp)? sin }(42— 8) 
exp(iA®)dA + exp{iky) | exp(zA”) (32) 
- 0 


oO 


which is similar to the solution given by Lamb (1942, p. 540). 
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6. Discussion 

The discussion will be concerned only with ¢,, which represents the 
rotational effects. The part of the solution ¢,+¢, is well known and has 
been fully investigated in the past. 

The expression in square brackets in equation (31) occurs in other 
diffraction problems, but the tabulation of it appears to be confined to a 
short table in the range 0 (-01) 1-0 for |x| and 0° (1°) 45° for arg s (Clemmow 
& Mumford 1952). The function would be determined for all values of 
arg z by a table covering the range —45° <args<45°. ‘To obtain a general 
idea of the behaviour of ¢, near the origin, further values of the function 
have been computed. In figure 4, the real and imaginary parts of ¢, are 


@=0 


/ 


6 +307° 


\ | | ° 
| J 
os 
3 


/ 


Figure 4. The real (dotted line) and imaginary (full line) parts of 
ls] o1 | for s=3/5 at different angles @ to the barrier. 


plotted against p for lines radiating from the barrier’s edge at different 
angles to the barrier for s=3/5 (p=4,5). ‘The amplitudes are given as the 
ratios of ¢, to |¢,]. ‘The horizontal scale is in units of (gh)'? y, which, for 
the latitude and depth of the North Sea, means units of 200km approxi- 
mately (y=10 ‘sec '). If the asymptotic behaviour of ¢, for large values of 
p is investigated, it is found that, except near the barrier for y 0, ¢; is of the 
order of 1/p!* for large p. Near the barrier, the asymptotic behaviour for 
v>0 is 


~ 2ia( exp(tkx cosh — ky sinh 8)+O( 


~ = exp{(tox yy) o( 3): (33) 


| 
by | \ = 
\ 
\ 2 'e wo | t2! 4 
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‘hus only in the region just behind the barrier are waves propagated which 
do not die out away from the edge of the barrier. It follows from equation 
(33) that, when p -3 5, these secondary waves are, at large distances down 
the barrier, larger in amplitude than the incident waves, although they die 
out exponentially away from the barrier. 

For waves to be propagated into the region behind ‘it, the barrier must 
lie in the right half of the v-axis, as in the problem studied. _ If it is in the 
left half (v--0), we may expect waves progressing down the front of the 
barrier. ‘his leads to the interesting possibility that, if the barrier is of 
finite length (an island in mid ocean), a certain amount of energy will be 
trapped in the form of a wave progressing round the barrier in a clockwise 
direction. 


I am indebted to Mrs P. Edwards for the computation of the integral 
in equation (31). 
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SUMMARY 

This paper considers the two-dimensional laminar boundary 
layer on an infinite flat plate normal to an oncoming stream, the 
plate making transverse oscillations in its own plane. The exact 
solution is shown to depend on a single ordinary differential 
equation, containing the frequency of oscillation as a parameter. 
Series methods are employed to evaluate the solution for small 
and large values of the frequency, and enough terms are calculated 
to give the solution with satisfactory accuracy over the whole ~ 
frequency range. It is observed that the solution satisfies the full 
Navier-Stokes equations. 

For a cylinder of arbitrary section making transverse or 
rotational oscillations, it is shown how the results for a flat plate 
can be used to describe the boundary layer in the neighbourhood 
of the front stagnation point. Difficulties which arise in extending 
the solution to cover the remainder of the cylinder are discussed, 
and an estimate is made of the fluctuating torque on a circular 
cylinder making transverse oscillations. 


1. INTRODUCTION 


Attention has recently been paid to the laminar boundary layer in 
flows in which there are fluctuations superposed on a basic steady low- 
speed motion. With the two-dimensional flow about an infinite cylinder 
as the basic motion, Wuest (1952) has considered the effect of oscillations of the 
cylinder in the axial direction, and Lighthill (1954) has considered oscillations 
parallel to the oncoming stream. However, in flutter problems, rota- 
tional or transverse oscillations are of chief importance, and it is with these 
cases that this paper is concerned. From the point of view of the boundary 
layer, it is immaterial whether it is the cylinder or the oncoming stream 
which is oscillating, as the inertial effects of an acceleration applied to the 
whole system are countered by a uniform pressure gradient, and the relative 
motion is completely unaffected. Compressibility effects will not modify 
this conclusion provided that the acoustic wave-length is large compared 
with the dimensions of the cylinder section. 

The results of this paper all follow from a detailed analysis of the two- 
dimensional flow against an infinite flat plate normal to the stream, the 
plate making transverse oscillations in its own plane. In §2 it is shown 
that for all amplitudes and frequencies of oscillation, the perturbation of 
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the flow is given by a single ordinary differential equation, containing the 
frequency as a parameter. Series expansions are obtained, valid for small 
and large frequencies, and enough terms are calculated to give the solution 
with satisfactory accuracy over the whole frequency range. Further, the 
solution satisfies not only the boundary layer equations, but also the full 
Navier-Stokes equations. 

At the lowest frequencies, the problem is identical with that of a steadily 
moving plate, and the perturbation to the boundary layer profile is deter- 
mined simply as the derivative of the basic profile, the well-known Hiemenz 
function. At the highest frequencies, the perturbation is a shear-layer, 
exactly as on a plate oscillating in a fluid at rest. It is found that the phase 
of the oscillating component of the skin-friction is always in advance of the 
velocity fluctuation, the phase-advance being proportional to the frequency 
for low frequencies, and approaching the limit }7 for high frequencies. 

This behaviour is similar to that found by Lighthill (1954) for oscil- 
lations in the magnitude of the oncoming stream. Although he had avail- 
able only the first term in each of the series for large and small frequencies, 
he was able to show that these are sufficient to cover the whole frequency 
range in that case. Stuart (1955) has considered the case of a fluctuating 
stream parallel to an infinite porous plate, and has obtained a solution 
which is valid for all frequencies. ‘The phase-advance again has the same 
qualitative behaviour, but, as in the present case, the first terms in the low 
and high frequency approximations are not themselves sufficient to cover 
the whole frequency range. In the case of axial oscillations of a cylinder 
considered by Wuest (1952), it is well-known that the axial motion has no 
effect on the flow in planes normal to the axis, and that the axial flow satisfies 
a linear differential equation. As well as considering more general cases, 
Wuest solved numerically the equation for the flow near the stagnation 
point for one particular frequency. It may be noted that, for general values 
of the frequency, the solution could be found by a straightforward appli- 
cation of the methods of the present paper. It is clear that the variation of 
phase-advance with frequency is similar to that in the other cases which 
have been treated. 

In §3,it is shown how the results of §2 may be applied to the boundary 
layer near the stagnation point on a cylinder making transverse or rota- 
tional oscillations. ‘l‘he first obstacle to be overcome in a consideration 
of the boundary layer further downstream on the cylinder is the difficulty 
of specifying the form of the fluctuating component of the velocity distribution 
at the edge of the boundary layer. In the case of streamwise oscillations, 
Lighthill was able to assume that the velocity distribution outside the 
boundary layer remained unchanged in form, and he succeeded in developing 
a very satisfactory Karman—Pohlhausen treatment of the boundary layer 
on a cylinder of arbitrary section. In the present case, it may perhaps be 
hoped that the considerations of § 3 will serve as a starting point for a Karman— 
Pohlhausen attack on the problem, and enable the solution to be continued 
away from the awkward region around the fluctuating stagnation point. 
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2. OSCILLATING PLATE 

The problem to be considered in this section is that of the two- 
dimensional flow of an incompressible fluid against an infinite plate normal 
to the stream. The theory will be developed for the case in which the 
plate makes harmonic oscillations in its own plane, and it will then be 
shown how the results can be immediately applied to the case where the 
plate is fixed but the point to which the dividing streamline in the oncoming 
flow is directed oscillates. 

Take Cartesian coordinates (x, y) fixed in space, the x-axis being along 
the plate and the y-axis normal to it, so that x=0 is the dividing streamline 
in the steady flow outside the boundary layer on the plate. Thus, if (w, v) 
are the corresponding components, outside the boundary layer u= U=cx, 
where c is a constant, and so, by continuity, v= = —c(y—5). The value 
of 6 will be determined in the course of the solution, and found to be constant. 
U and V are seen to be the velocity components in the two-dimensional 
potential flow against the plane y=6. ‘The plate is assumed to have velocity 
ae in the x-direction, where a and w are constants. 

The boundary layer equations are 


at “ax dx OP (1) 
and 
Ou dv 
where v is the kinematic viscosity. These are to be solved with the boundary 
conditions 
u=ae™, v=0, aty=0, and u— U=cx as y—> o. (3) 


For a=0, the plate is fixed and the solution takes the classical form 
obtained by Hiemenz in which the velocity components are u=cx/f'(y), 
v= —(cv)!?f(y), where 7 = (c v)!?y,a dash denotes ditferentiation with respect 


to n, and f satisfies the equation 
+ff' -f? +1=0 (4) 


with boundary conditions f(0)=0, f’(«)=1. Equation (4) has 
been studied in detail by many writers. A ten-figure tabulation of f, f’ and 
f”, performed by Dr N. E. Hoskin on the Manchester electronic computer, 
proved invaluable during the subsequent calculations of this paper. 

To satisfy equation (3), we look for a solution in the form 


u=cx f'(n) + aed(n), v= —(cv) f(y), (5) 


where the meanings of f and 7 are unchanged. Since 7 is a function of y 
only, the continuity equation (2) is still satisfied. Equation (1) now becomes 


G2 
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The first bracket is zero, by equation (4), so both the equation and the 
boundary conditions (3) will be satisfied provided 


Thus, as for a fixed plate, there is only one ordinary differential equation 
to be solved. tis to be noted that this result holds for all values of w and a, 
so there are no restrictions on the amplitude or frequency of the oscillation. 
Further, this solution, like the Hiemenz solution itself, does in fact satisfy 
the full Navier-Stokes equations and not just the boundary layer equations. 
The only term omitted in the momentum equation (1) is 0?u/dx?, which 
vanishes, and, since there is no extra contribution to v, the momentum 
equation in the y-direction is unaffected. 

Since f(7) ~ 7-0-6479 as »-» x, the potential flow outside the 
boundary layer is that against the plane y = 0-6479(v/c)#?. 

The only parameter in equation (7) is the frequency ratio w/c. ‘Two 
series solutions will be developed, valid respectively for small and large 
values of wc, and it will be shown that sufficient terms are here calculated 
to enable the whole frequency range to be covered satisfactorily, so far as 
the value of the skin-friction is concerned. 


Small values of w/c. 
Consider first the case w =0, which implies that the plate velocity has the 
constant value a. ‘hen, by equation (7), 6=¢», where 


Fo(0)=0. (8) 
Now, differentiation of equation (4) gives 
(9) 


and so ¢y=/f" satisfies equation (8). Also f”(0) = 4 = 1-2326, a value found 
by the integration of equation (4), and f’(«)=0. Hence 


1 
by= Gf” = (10) 


satisfies both the boundary conditions. ‘Thus, in the flow against a plate 
moving with steady velocity a, the velocity components have the simple form 


(11) 
For small but non-zero values of wc, write 
20 n 
> dun. (12) 


Consideration of the terms in in equation (7) shows that, for n> 1, 


bt $n(0)=0, (13) 


5 
| 
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Equation (13), like equation (8), has f” as one integral belonging to its 
complementary function. The other integral is easily found to be 


f "| a dy, where f “=| fd. The method of variation of parameters then 
0 0 


leads to the solution 


fe bys do (14) 
/0 
where 
ms 
“Le 


Further details of this solution are given in the appendix. It is shown that 
the boundary conditions are satisfied, and that ¢,, tends to zero exponentially 
as 7 tends to infinity, for all n. 

Actually, ¢, can be obtained directly in a rather simpler manner. ‘The 
equation for 4, is 


6:00) =0, (15) 


A particular integral is seen to be ¢, = af To this must be added suitable 


multiples of the complementary functions f” and f”/ so as to satisfy the 
boundary conditions. Equation (56) of the appendix shows that the 
required solution is 

= (16) 


‘This expression can also be obtained from equation (14), with the help of 
equation (9) and its derivative. 

Numerically, we shall confine our attention to the value of the skin- 
friction, equal to u(du/dy),-9. The oscillating component of the skin- 
friction 7 is given, from equation (5), by 


T cv\i2 
inth'(()). 
(17) 
The contributions to ¢'(0) from 4) and ¢, are immediately available, since 
it is known that f’(0) =A, f”’(0) = —1, while that from ¢, require only the 


computation of | fed, dy. Numerical integration yielded the value 
0 
—0-1167, so that the first three terms of the series give 


2 
=0-8113 + 0-4932°" +0-09472 (18) 
| c 


(it is natural to consider — 4’(0), since a negative value of the skin-friction 
would be expected to be associated with a positive velocity of the plate.) 
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Large values of w/c. 
When w ¢ is large, it is convenient to change the variable in equation (7) 
from 7 to 


Y = (tw, = (tw, (19) 
d d 
Write (c/tw)'?=2 ; then, since equation (7) becomes 
p**—d=al fed —fd*} : d(0)=1, d( 0)=0, (20) 


where a star denotes a derivative with respect to Y. 
The expansion for near 7 ts 


f= 5 An 6” 120 A n 


2 6 2520 
For large frequency, ~ is small, so we look for a solution in series by writing 
d= (22) 
n=0 


and equating to zero the coefficients of successive powers of x in equation (20). 
The boundary conditions are 


w(0)=1, 9,(0)=0 (n>1), (all n). (23) 
The equation for @p is 
(24) 
and the required solution is 
(25) 
This is the familiar shear-wave solution due to Stokes, which gives the flow 
due to a plate oscillating in its own plane in a fluid at rest. 
‘The second and third equations show that q, and g, are zero. The next 
two equations are 


Now, when > @,¥*, 
k=0 
so the solutions of (26) and (27) satisfying the boundary conditions (23) are 
and 
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The function q; is identically zero. In the equations for ¢,, y; and qx, 
there are contributions to the right-hand side from gy; and q, as well as from 
Go. Making use of equation (28), we finally obtain 


= — Ae {By + y24 y34 ys 
+ "+ (32) 
Since 


ao 


$'(0)= =4e(0) = 2” 


the contribution to the skin-friction from the calculated terms is given by 


3_ + 


8 16° 128 512” 

1 fe\-¥2 3y2fc\32 33, 139\/2/c\*?) 

129 139\/2/c\72 
256 - (<) 


since 4 = = (1 


‘The variation of the real and imaginary parts of —4'(0) with frequency, 
as given by (34) for high frequencies and (18) for low frequencies, is shown 
in figure 1. The close agreement between the two formulae in the region 

“near w c= 1 in each case affords good evidence that the terms calculated are 
sufficient to predict a reliable value of the skin-friction over the whole 
frequency range. The fact that the real and imaginary parts of —¢’(0) 
have the same sign indicates that, in all cases, the skin-friction is advanced 
in phase relative to the velocity fluctuation. The variation of the phase- 
advance angle @,, which is the argument of —¢'(0), with frequency is also 
shown in figure 1. 4, rises steadily from 0 to z 4, as the frequency increases. 
For small frequencies, #,=0-6078 c+ O(w* indicating a time of 
anticipation #, w =()-6078 c, which is independent of frequency as long as 
cis negligible. 

When the plate velocity is not sinusoidal, but may be expressed as a sum 

of sinusoidal components, it is clear from the form of equations (5) and (6) 

that the combined effect is the sum of the effects of the individual components. 


pepe 
4 
- 
bis: 
a 


104 M. B. Glauert 


In particular, if the frequencies of all the components lie in the range for 
which w?/c? is negligible, the time of anticipation in the fluctuating com- 
ponent of the skin-friction remains as 0-6078/c. 


j 37083 
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Figure 1. Oscillating plate. Variation with frequency of the real part (#) and the 
imaginary part of and of the phase-advance angle 44. 


Oscillating stream. 

When the dividing streamline of the oncoming stream oscillates in 
position, but the plate is at rest, the situation differs from that already con- 
sidered only by the superposition of a uniform, though not constant, trans- 
verse velocity, which has no effect on the relative motion. It follows that 
the solution obtained above can be applied at once to this new case, though 
the details require a little care. 

With coordinates as before, the velocity just outside the boundary layer is 


(35) 


The stagnation streamline is thus instantaneously directed towards 
x= —(b/c)e. The plate itself is at rest. 

Take a new origin of coordinates O, at x= de", and write x, =x — de as 
the new coordinate along the plate. ‘lhe point O, has velocity twde*”, so the 
stream velocity at x, relative to the new axes, is 

U,=cx + bel — inde = ex, + (b + cd—iwd)e™*, (36) 


Now d must be chosen so that l', does not fluctuate ; hence 


h+d(c—iw)=0. (37) 


\ 
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In the new axes, the plate velocity is —iwde“, and on writing 
a= —iwd (38) 
we recover the case of the oscillating plate precisely as treated above. 
The velocity distribution is 
= cx, + (39) 
On expressing x, and a in terms of x and 4, we obtain the velocity in the 
original coordinates as 


=U, +twde™ = cx f'(n) + be y(n), i (40) 
where (41) 


Alternatively, it may be easily shown directly from the boundary layer 
equations that x satisfies 


Equation (42) is more troublesome to consider directly than was equation (7). 
But it is easy to check that the value of y given by (41) does satisfy equation 
(42) for all values of w c. 

The result (42) is again exact, with no restrictions as to amplitude or 
frequency. 

For small values of w c, since 

c 

equation (41) gives 


xa f+ SU — 14 byt tbs) +o 
(43) 
and for large values of w c, since 
249, +..., 


As in (17), the oscillating component of the skin-friction is given by 


\ 


= ( (45) 


The variations with frequency of the real and imaginary parts of x’(0), and 
of the phase-advance angle @,, are illustrated in figure 2. For small values 


E 


(67) 


of the frequency, a,=(1 ~ -a)2=03482, indicating again a constant 


time of anticipation 6,,w=0-3418c. At higher frequencies, 4, rises 
steadily towards the limiting value }z. 


vin 
| 
Pig 
Pe 
= 


106 M. B. Glauert 


Thus, in each of the cases considered here, the behaviour is qualitatively 
similar to that found by Lighthill for the case in which the stream velocity 
fluctuates in magnitude. For the Hiemenz layer, which corresponds to 
the plate making oscillations in the y-direction, the time of anticipation at 
low frequency was determined as 0-18 c, considerably less than in either of 
the present cases. Lighthill’s analysis was restricted to infinitesimal 
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Figure 2. Oscillating stream. Variation with frequency ot the real part (P) and the 
imaginary part (7) of y(0), and of the phase-advance angle 6 


disturbances of the steady stream. ‘This was unavoidable, since c fluctuated 
during the motion, and hence the additional velocity had components in 
both the v and y directions. It may be noted that he found it possible to 
join up satisfactorily the high and low frequency formulae, although he had 
available only the first fluctuating term of the series in each case. 

Finally, a word may be said on the question of fluctuations in heat 
transfer from a heated plate. This problem was discussed in detail by 
Lighthill for his case. For the Hiemenz layer, the temperature 7 is a 
function of y only, satisfying 

of 


T=T, at} [--T, as4 x (46) 


where « is the thermal diffusivity, and 7), and 7 are the constant tempera- 
tures of the plate and stream. Since, for the problems of this section, the 
oscillations do not involve any change in 7 or in the boundary conditions, 
the steady solution of equation (+6) continues to apply unchanged. In 
other words, the movement of the plate has no effect at all on the tempera- 
ture distribution and the rate of heat-transfer. 
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3. OSCILLATING CYLINDER 
In this section, we consider how the results of §2 may be applied to a 
cylinder of arbitrary cross-section. 
Any two-dimensional oscillatory motion can be considered as a com- 
bination of the following five basic motions: 


(i) cylinder fixed, stream oscillates in magnitude ; 

(ii) cylinder fixed, stream oscillates in direction ; 
(iii) stream constant, cylinder oscillates in the stream direction ; 
(iv) stream constant, cylinder oscillates in the transverse direction ; 
(v) stream constant, cylinder oscillates about its axis. 


Here, the word ‘stream’ is used to denote the fluid velocity at large 
distances from the cylinder. 

Cases (i) and (ii), which are equivalent, were analysed in detail by 
Lighthill. He assumed that the velocity U(«) outside the boundary layer 
remains unchanged in form, merely changing in scale in proportion to the 
velocity of the cylinder relative to the stream. 

In the general case, the flow outside the boundary layer differs in form 
from the corresponding steady flow, except possibly for a circular cylinder, 
and cannot be determined by boundary layer theory alone. A first approxi- 
mation, valid for low frequencies, would be to assume the external flow to 
be the quasi-steady flow, i.e. that due to a steady stream given by the in- 
stantaneous relative velocity of the cylinder and the distant fluid. Cor- 
rections due to discrepancies in the separated region at the rear of the 
cylinder might be taken into account in a second approximation. If the 
cylinder is of aerofoil section, with circulation, the question of the fluctuating 
component of circulation provides another problem. Except at small 
values of w c, vorticity shed at the rear will still be in the neighbourhood of 
the cylinder after a few periods, and so will have its effect on the velocity 
distribution. For, if |” is a representative velocity and / a representative 
length of the cylinder, so that |’ = cl, it is seen that the distance travelled in 
a period is 2z/c'w. For small values of w ¢, the quasi-steady approximation 
* would still be appropriate. For large values of w c, the best approximation 
would be to assume that the circulation remains constant, even though the 
Kutta—Joukowski condition is no longer satisfied. ‘These considerations 
indicate that Lighthill’s treatment is inadequate for an aerofoil with cir- 
culation, even if the fluctuation is confined to the magnitude of the on- 
coming stream. 

However, it is always true that the velocity over the surface near the 
stagnation point can be written as LU =c(x—wxy), where c and x, have 
fluctuating components. If these fluctuations are small, the contributions 
due to the variations of c and x, may be added separately to the basic steady 
flow. The effect of variations in c are covered by Lighthill’s analysis, and, 
with c constant, the analysis of §2 becomes applicable. As discussed above, 
the value of x, may depend on other than purely quasi-steady considerations, 
We now consider the cases (ii), (iv) and (v) in turn. 
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Case (ii) calls for an immediate application of the oscillating stream 
formulae of $2. If the stagnation point is at x=se', x being measured 
round the cylinder perimeter, and if c is the velocity gradient at x= 0 in the 
undisturbed flow, then, by equation (35), the solution is given by substi- 
tuting b= —cs in equation (40). ‘This holds whether or not the value of s is 
that given by the quasi-steady flow at the appropriate incidence. It is 
perhaps best to consider s as being found by experiment. The extra 
contribution due to any fluctuation in ¢ is given by Lighthill’s formulae. 
It is to be noted that, to the first order in the amplitude, no such contribution 
arises in the case of a symmetrical cylinder when the flow oscillates about 
the axis of symmetry. It is likely that, for nearly circular cylinders, the 
contribution to the oscillating skin-friction deduced from equation (45) is 
indeed the major part for the whole cylinder, since over the shoulders of the 
cylinder, where all resemblance to a linear gradient has gone, the fluctuations 
in the quasi-steady external flow are much reduced. In any case, the 
boundary layer is likely to have become turbulent by this stage. 

Case (iv) might, at a first glance, be thought to require the application of 
the oscillating plate analysis of §2, but this is not so. ‘The essential point 
is that, whereas there the axes were fixed in space, here it is natural to consider 
them as fixed in the cylinder. ‘This is equivalent to superposing a uniform 
though fluctuating velocity on the system, which has no effect on the relative 
motion, and so case (iv) is immediately reduced to case (11), with an oncoming 
stream as given in the new coordinate system. ‘his procedure is identical 
with that by which case (iii) is related to case (1). 

Case (v) may also be reduced to case (ii), by choosing axes fixed in the 
cylinder. According to the boundary layer approximation, the only effect 
of rotation, as of curvature, is a small but still negligible pressure change 
across the boundary layer. ‘The flow L(x) at the edge of the boundary 
layer will have a different value from that in the corresponding situation in 
case (11), but may be found in suitable cases by quasi-steady considerations. 
For a circular cylinder, an alternative approach is to use fixed axes and to 
apply directly the oscillating plate analysis of § 2, the notation being suitable 
as it stands; a more detailed study of this case, with particular attention to 
a steadily rotating cylinder, will be made in a later paper. 

As a final example, we attempt a quantitative estimate of the overall 
effects on a circular cylinder making small transverse oscillations. We 
assume that the flow at the edge of the boundary layer is quasi-steady, and 
can be written in the simple form 

U=c(x+%), x, (47) 
where x=, is the stagnation point and x, is a constant. We assume that 
the fluctuation effects are negligible outside this region. ‘Typical values 
for a cylinder of diameter d in a stream IV are c=3-6V d, x, =0-4d. If the 
transverse velocity component is Se", then, on the quasi-steady hypothesis, 

Pd 
w= 


tor small values of 3, 
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For |x+x9| < x, the velocity in the boundary layer is given by equation 
(40) as 
u=cxf' + cxox 

= (x + Xo) + (48) 
The first term of (48) is a contribution to the quasi-steady flow round the 
cylinder, which results in a fluctuating lift force of amplitude BD V’, where 
D is the drag force on the cylinder, in phase with the velocity fluctuation. 
The second term of (48) involves a skin-friction given by 


T 
and the values of this expression for all values of w are given in $2. Now 
x'(0)— 4 =0-4213 + (=): (50) 


and thus, for small values of w c, the skin-friction acting over the region 
|x +x] < x, gives rise to a fluctuating torque about the cylinder axis, in 
phase with the acceleration, of amplitude per unit span 7 given by 


T =0-4213 pBewx,d(v c)"2. (51) 


For the typical values of ¢ and x, given dbove, and for the particular values 
w=0-5¢e, B=0-1V, (51) becomes 


» \ V2 
T= 0-016( 73) (52) 


Since the fluctuations in velocity for |x +x | > x,, assumed zero here, will 
in fact be much smaller than those for |x + x)| < x,, it seems likely that (51) 
is indeed a reasonable estimate of the torque on the cylinder due to skin- 
friction. 

For the largest values of w c, y'(0)— A ~ (iw,c)'?, but before this regime 
is reached the basic theory will in many practical cases have broken down, 
as the acoustic wavelength will have become comparable with the cylinder 
diameter. 


APPENDIX 
Equation (13) of § 2 was ‘ 


+fb,' —f'bn=bn-13 $,(0)=0, ¢,(0)=0 (n>1). (53) 
The integrals of the complementary function are f’and f”/, and it is clear 


that e diverges as x. ‘To investigate convergence and 
0 

to obtain numerical coefficients, it will be necessary to consider the 


behaviour of J in some detail. 
The integration of equation (+) has shown that as 7 -- %, 
f ~ €=n-0-6479, ~ +0-1496, 


Hence et ~ Be 


where B=0-8610. 
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Writing f=&+g in equation (4), and ignoring terms in g?, we obtain 
g. q g 
— 20 (). (54) 
Two differentations give g’+ég'¥=(0, and hence for 
some constant (. Integrating, we have 
e~™ dx, (55) 


the well-known error function of which tables are available. Now g’”=f"", 
and use of the tabulated values of f, together with equation (4) itself, gives 
C=0-6451. 


As > ~ Ce*™, 
Hence [~ Stel dé 
é (56) 
and 


Returning to equation (+7), where d,_, is supposed known, we obtain a 
solution of the form 


o,=h, f’ +k, f'l, (57) 


provided that 
h' ,, ki 


a 


In view of (56), the boundary conditions on ¢,, require 


n 


| h, dn, k,=—| k, dy. (58) 
“0 
Thus 4, has the form already given in equation (14). 

It remains to be verified that 4, does in fact tend to zero at infinity ; 
indeed from general boundary layer experience we would expect ¢,, to tend 
to zero exponentially. ‘his is the case, as is shown by the following 
induction argument, 

Suppose that for some ¢,.,=O (e-*') for large where 
€ =7 —0-6479 as before, and here and in what follows O (e ?*’) is taken to 
include O(ge ”*). ‘Then, as a consequence of the results obtained 
above, h’,=O(1), ki =O(e **), and so by (57) and (58), 6, =O(e-*"). 
Since dy =O(e~*"), the result is true for all by induction. 


A 
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SUMMARY 

The region of finite vorticity near the edge of a diffracting 
wedge is investigated. Dimensional analysis gives the dependence 
of the circulation and the velocity of the vortex region on the pulse 
strength. A close estimate of the magnitude of these quantities is 
obtained by replacing the vortex region by a single concentrated 
vortex. ‘The theoretical conditions at the sharp edge are dis- 
cussed and compared with observations of real fluid behaviour. 
A short account of the theory of the core cf the spiral vortex 
sheet in a perfect fluid is appended. 

‘The research reported herein was supported by the United 
States Air Force, through the Office of Scientific Research of the 
Air Research and Development Command. 


1. FORMULATION OF THE PROBLEM 

Observation of the diffraction of shock waves by sharp wedges shows 
the existence of a region of finite vorticity near the edge. ‘The pheno- 
menon is shown in plate 1, which reproduces one of several shadowgrams 
taken at the University of Michigan (Uhlenbeck 1950). The purpose of 
the present paper is to estimate the position, the velocity and the total 
circulation of the vortex region, in the limiting case of weak shocks. 

On using the acoustic approximation for weak pulses, the classical 
diffraction problem of Sommerfeld is obtained if the generation of vorticity 
is ignored. The solution of this problem for non-zero wedge angles was 
given by Friedlander (1946) and, in a modern and strongly simplified form, 
by Keller & Blenk (1951) and Miles (1952). ‘This solution exhibits (in 
general) an infinite velocity near the sharp edge; this fact indicates why the 
vortex-free solution is not found in reality. 

Experience shows that instead of flow past the sharp edges with infinite 
velocity, a discontinuity sheet forms, originating at the edge, and the velocity 
at the corner remains finite. Viscosity must be essential for this deviation 
from the vortex-conservation laws of inviscid flow; however, the problem 
seems amenable to the now-classical aerodynamic technique in which the 
detailed mechanism of vortex production is by-passed and the flow with 
vortex sheet treated by inviscid fluid theory. To do this, the Kutta 
Joukowski condition is formulated: vorticity is produced at any instant at 
such a rate that the resultant flow, with the velocities induced by the vortex 
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elements taken into full account, possesses a finite velocity at the corner. 
With this condition a problem is obtained which, at least for incompressible 
flow, is soluble in principle, although the actual calculations may lead to 
formidable difficulties. (Arguments about the existence of such solutions 
will be given later.) Qualitatively, however, by considering the motion of 
newly produced vorticity under the influence of the vorticity generated 
previously, it is easily understood that the discontinuity sheet will be wound 
into a narrow spiral, outside of which the flow may be taken as irrotational. 
Compressibility effects will complicate the phenomenon further, especially 
towards the core of the vortex; this point will be discussed later. 

On the basis of the foregoing arguments, itis proposed to treat the problem 
by neglecting viscosity (and heat conduction). ‘This assumption has strong 
support from experimental evidence. In inviscid flow, the diffraction 
pattern of a sharp pulse (of any strength) at an infinite wedge is exactly 
similar to itself for all times t. With a coordinate system x, y centred at 
the vertex, the velocity, the pressure, and other functions of state become 
functions of x tf and yt only. ‘This ‘ quasi-steady’ or ‘ conical’ property 
of the flow field is strikingly well realized, as flow pictures taken at con- 
secutive times show, for diffraction patterns with vortices present. Scrutiny 
ot the experimental data reveals only a small deviation from quasi-steady 
flow, which is accounted for by viscous effects, and for the purpose of the ~ 
present theory the inviscid assumption appears well justified. 

Next, in view of the restriction to weak pulses, it is proposed to treat the 
diffraction problem with vortex separation in the linearized approximation, 
i.e. a potential 4 will be assumed, satisfying the wave equation 

2, 
(1) 


ox? 


(a=velocity of sound, a constant). For quasi-steady solutions, any 
derivative of the potential 4 will be a function of =x at and »=y at 
only; therefore, ¢ has the form 


b=at f(§,). (2) 
Introduction of these variables into (1) leads to the equation 
fl-&)+f,,(1 — 1°) — 2&nf., =9. (3) 


Observations show (see plate 1), and the subsequent calculations 
verify, that, for weak shocks, the vortex region is limited to a domain whose 
radius is small compared with the radius at of the diffraction region. Here, 
|€|<1 and |x|<1, so that (3) may be well approximated by Laplace’s 
equation, 


fest = 9. (4) 

Thus the acoustic approximation leads to incompressible flow in the 
central region of a quasi-steady flow field. 

The incompressible treatment of the central region appears well justified 

for the problem with vortex generation, as the resultant velocities remain 

finite due to the superimposed effect of the vortex sheet. Experience 


Me, 
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Plate 1. Shadowgram of a shock wave diffracted by a 90° wedge. M,—0-27, 


t—152 microseconds. Ahead of the shock, the air is under atmospheric 
conditions. (Reproduced, by kind permission of Prof. O. Laporte, from a 
report prepared at the University of Michigan.) 
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shows, however, that appreciable density changes occur towards the core 
of the vortex spiral. ‘These density changes were investigated experi- 
mentally and theoretically by Howard & Matthews (1955). From their 
results it may be seen that, in the limit of weak shocks, the density variation 
vanishes sufficiently rapidly to justify an incompressible treatment. 
Furthermore, for moderate shock strengths, density changes near the vortex 
do not necessarily invalidate the subsequent estimates based on incom- 
pressible flow. In order to be able to calculate the vortex motion under 
the influence of self-induced velocities in the presence of the wedge by 
means of the approximation of incompressibility, we need assume only that 
the Mach number of the flow velocities along the wedge surface remain 
small compared with unity. 


2. SOLUTION BY DIMENSIONAL ANALYSIS 
‘The linearized potential ¢ for the diffraction problem with vortex 
separation, as a solution of the wave equation (1), can be considered as a 
sum of two solutions : 


bp, (5) 


where d, is the vortex-free (Sommerfeld) solution, and <J,. is an additional 
part due to the existence of the vortex sheet. ‘The potential 4, is known; 
it depends on the strength of the pulse, the incidence of the oncoming wave, 
and the wedge angle. ‘Che dependence on the strength of the pulse can be 
given immediately. Let u, be the gas velocity in the rear of the undisturbed 
incident pulse. (‘This ‘ afterflow velocity’ is taken positive in the direction 
of the incident wave propagation.) All velocity components derived from 
the solution will be proportional to u,, so that one can write, with a simple 
renormalization of (2), 


ps = t f(E, 1), (2a) 


where f, is dimensionless and depends on the flow geometry, i.e. on the 
wedge angle and the angle of incidence enly. 

The additional solution 4, naturally also depends on uy and the flow 
geometry ; however, its dependence on uy, will not be as simple as that 
expressed by (2a). It is clear that the position of the vortex relative to the 
‘Mach circle’ (with radius at) will depend on the pulse intensity m7, and 
thus the velocities derived from ¢, will not be simply proportional to up. 
Using the fact that for u,<a the vortex region is very near to the vertex of 
the wedge, the dependence of 4, on uy will be found by dimensional analysis. 

The function f, 1n (2a) is a solution of equation (3), and is known. For 
|y|<1, must become a solution of Laplace’s equation (4). 
Therefore, 

PF ALF SMO), (6) 


for |¢|<1. It is easy to write down the form of the function F(¢), as it 
must represent the complex potential for incompressible flow past an infinite 
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wedge. Let the wedge with the angle £ be formed by the positive real axis 
and ray with the argument —f. ‘The function F, takes the form 
(7) 


where 


(8) 


so that, for 0<fB<7, m lies in the range }<n<1. The dimensionless 
quantities A,, A,. .. become real for the wedge position specified above ; 
they depend on the flow geometry and can be given only if the Sommerfeld 
problem is solved, i.e. if /, is found as a solution of the full equation (3). 
Expansion of the Sommerfeld solution around ¢=0 naturally confirms the 
form of (7), but only the first two terms represent analytic functions of the 
present variable ¢. This means that if more terms were needed for the 
treatment of the vortex generation problem, the use of the approximate 
equation (4) would be unjustified. 

As | <n—1, only the first term in (7) represents a ‘singular’ flow past 
the wedge with infinite velocities, which is taken as the ‘cause’ of vortex 
separation. It may be imagined that the whole phenomenon of vortex 
generation is not only ‘caused’, but also fully ‘dominated’ by the first 
term of (7), which would then be the only one that need be kept for the 
subsequent analysis. This assumption will be made now, and its justifi- 
cation will be sought a posteriort. 

Returning to dimensional variables, put 


=auy t Fy, z=atl, (9) 


~ n ~ \2n 
(=) “\at J 


A, n at A, 2n 


so that 


(10) 


If only the first term plays a role, then there is only one physical constant in 
the problem: A,, which has the dimensions of a velocity to the power 2—n. 
Now consider the vortex region, which is ‘generated and dominated’ 
by the first term of the flow potential (10), and let it be characterized by the 
total circulation I’, the position of the core (centre of the spiral) z,, and the 
velocity w, of the point z,. As the resultant flow is also quasi-steady, the 
time-dependence of these quantities is given as follows: 


Sr Wp x 0°, ee. (11) 


The dimensional factors multiplying these time-dependencies can be formed 
only by the quantity A,, and in only one way: 


{ 
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or, using the expression for A, in terms of a, uv, and A,, 


u 1(2-7) ) 
zp oc — at, | 
a 
1(2-n) 

Up | 

Wp Ale a, (13) 
— = — aug t. | 
a a 


‘he use of the factor 4} " does not make much sense at present, as a 
further factor dependent on the flow geometry is implied in (13). But (13) 
gives the correct dependence on the ‘afterflow Mach number’ uw, a (- M,). 
which characterizes the pulse strength. Indeed, if @ and wu, are not inde- 
pendent constants of the problem but occur only in the combination K,, 
then (13) represents the only possible solution. 

The explicit dependence of 4, on VW, is given as d,, < I, so that 


b, Mi (2 au 


for |C|<1. 

The validity of the original assumption that the first term of the 
expansion (10) is predominant, can now be shown in an ‘iterative’ way. 
Introduce z;. from (13) into (10) and compare the magnitudes of the different 
terms. ‘The Nth term becomes 


Nn 
be (Zp) =autA,( 
bs wh r) 0 (2 


tA, ANC (14) 
It is seen that the magnitude of the terms diminishes rapidly with increasing 
\ if M, is small. Thus, in the limiting case of weak pulses, the 
vortex region will really be ‘generated and dominated’ by the first term. 
Howard & Matthews (1955) found that the corresponding dependence 
z, at was excellently verified by their own experiments. 

It must be emphasized that the simple result obtained for diffraction is 
not always valid for vortex generation problems; indeed, it is quite excep- 
tional. In wing theory there is considerable interest in problems of vortex 
separation from sharp edges. ‘lhe method generally used is the replace- 
ment of the vortex region by a single concentrated vortex, as will be dis- 
cussed shortly. Only one case in wing theory corresponds to the diffraction 
problem : for a rectangular thin wing in supersonic flow, the vortex separation 
from the side edge is perfectly analogous to the diffraction by a wedge with 
3=0. This case was treated by Cheng (1955), using the single-vortex 
model, and his results agree with (13). Other wing problems can be related 
to the diffraction case only if the wedge is supposed to move, with some con- 
stant velocity zy. This causes additional complications, since the pheno- 
menon depends strongly on the direction of the motion of the wedge. 
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‘The simplest cases are obtained for wedge motion in the direction of the 
bisectrix of the wedge angle, as otherwise the wedge motion causes vortex 
separation even without the diffracted pulse. As an example, the case of 
zero wedge angle (8=0, n=) will be considered. ‘The potential of the 
constant flow relative to the edge is simply zyz. If 7) < a, a development ot 
d. relative to the moving edge can be expected to have the same form, (10), as 
for 7, =0, except that the 4,, 4,, .. become functions of 7, uv), and there 
will be an additional term of the form 

Ad,=t9 2 
due to the proper motion of the plate. Now, the ‘iterative’ justification 
of the predominance of the first term in (10) fails, because of this additional 
term. From (13), with n=}, it follows that 


a 
Equation (14) yields, for 1, 
13 
a 


which is larger than the previous term only if 


A, u)\a 


If v is kept constant, this condition clearly fails in the limit u, - 0- 
Physically, this means that for a weak shock, and 7, > u,, the vortex is swept 
away from the edge by the velocity v7). If, however, 7, and uv, have the same 
order of magnitude, and wu, ay is still small, (13) becomes valid again. 

It may be noted that, for a wedge flow with ‘generalized’ quasi- 


oy, = (=) 


near s=(, a result corresponding to (12) may be obtained, if A, 1s the only 


steadiness such that 


important physical constant : 


L (16) 


3. ‘THE * SINGLE VORTEX’ APPROXIMATION 
In erder to find the numerical factors still missing in (13), use will now 
be made of an approximation mentioned previously: the whole vorticity 
region will be replaced by a single concentrated vortex. The results 
already obtained by dimensional analysis will incidentally be re-established. 
‘The complex potential near the vertex will be decomposed again into a 
vortex-free part and an additional solution due to the vortex: 


h=dbgt dp. (3a) 
According to the results of §2, put 
by = (10a) 


23 
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For the single-vortex model, ¢,. is the potential for a vortex of strength J” 
situated at the point z,=at@,., in the presence of the wedge, which must be 
a streamline of the ow(.4‘'dé!=constant). The solution is 


- 
on 
For this potential, the velocity vanishes as |€|~- ~<. It is easy to realize 


that this is the proper boundary condition at infinity. ‘The quantities ¢, 
and 4, are, respectively, the approximations for |¢|<1 to the solutions ¢, 
and 4, of the full wave equation. The Sommerfeld solution 4, takes care 
of all the boundary conditions along the Mach circle, imposed by con- 
ditions of proper joining to the pulses outside the Mach circle, these joining 
conditions are naturally the same for the solution ¢ (with vorticity) and for 
dy, so that the solution ¢,,=¢—d, fulfills ‘zero disturbance’ conditions 
along the Mach circle. ‘Thus, the solution 4, expresses the sole effect of 
the growing vortex region, and the vanishing of the solution on the Mach 
circle simply expresses the fact that this disturbance, originating at s=0, 
t=, can spread only with the velocity of sound. Evidently, as long as 
the incompressible approximation can be made, i.e. as long as z,< at, it 
makes no difference whether the ‘no-flow’ condition is fulfilled on the 
Mach circle or at infinity. 
From the total potential 
b= aly 4,(=) + (18) 

the complex velocity w, ‘at’ the singular point 3, will be derived. This 
is defined, according to Kirchhoff, as the regular part of the complex velo- 
citv at 2,, 


. 0d 
w,= lim as Sp. (19) 
‘The calculation of the limit, which requires some care, leads to the following 
result: 
any il’ fl—n cr 
We= A, — + }>. 20 
If the vortex is free, it moves with the Kirchhoff velocity w,, that is 
dz p 
=aCp. (21) 
dt ; 


(Note that the complex velocity is the conjugate of the physical velocity 
vector.) Equations (20) and (21) together give an important connection 
between J’ and Cp. 

A further complication arises from the fact that this single vortex with 
variable intensity is not free; strictly speaking, it is incompatible with the 
vortex conservation laws. The original vortex spiral naturally is free, i.e. 
it does not sustain any forces. By concentration of the spiral into a single 
vortex, its ‘umbilical chord’ has been cut, so that either the vortex strength 
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cannot grow, or, if /’ changes, an unbalanced pressure jump will be found 
on any line connecting the vertex of the wedge and the vortex. ‘This pressure 
jump is expressed by the fact that dé dt has a multivalued (logarithmic) term 
for variable /"; the magnitude of the jump is pd/* dt =constant. Due to the 
‘hydrostatic’ character of this jump pressure, the resultant force is the 
same on any line connecting the vertex and the vortex, and is equal to 
t; its direction is perpendicular to z,. 

It would seem that the single-vortex model has to be abandoned, due to 
its inconsistency with vortex-conservation laws. However, an ingenious 
proposal which enables this useful approximation to be retained has been 
made by Brown & Michael (1954) and by Edwards (1954): let the single 
vortex not be free, but subject to a Joukowski force, which cancels the un- 
balanced force on the vertex-vortex line. (Only an unbalanced moment 
will remain.) In order to present this condition in full generality, let the 
sheet be generated at a point zy, so that the (complex) force to be balanced 


1S 
The Joukowski force is 
dsp 
F, io( (23) 


[| dt dt 
It has been pointed out by Cheng (1955) that this equation may be used for 
any flows with vortex generation, even without similarity. In the general 
similar case expressed by (16),s,% ?” and x f?"!, so that (24) yields 
(with z,=0) 


or, for quasi-steady flow (m= 1), 


Wp = 2al p, (26) 


which replaces (21). It is seen that (21) and (26) ditfer only by a numerical! 
factor 2; the exact number to be taken is naturally unknown as long as the 
true spiral structure is undetermined, but (26) may be supposed to lead to 
a close estiinate. ‘The ‘free-vortex’ condition 


0; = = (27) 
at 


leads to the same result as (25) for m=}, when the single-vortex model 
becomes exact: [oc f?”" '=constant. ‘This interesting limiting case, which 
may be realized in special cases of incompressible flow, is certainly worthy 
of further investigation. 


‘ so that the condition ot the vanishing total force leads to 
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Returning to our problem, the elimination of w, from (20) and (26) leads 
to the following connection between ¢,- and I’: 


= Cp Cp = My Ay + (28) 
Putting ¢,. =o exp (#9), the imaginary and real parts of (28) yield, respectively, 
O0= M, A, no" sin nt — ——, 29 
al 
and M, A, cosnd — = cot mi. 30 
n 
Elimination of J” in (30) by (29) leads to the equation 
? 
- = M, A, o"(1—2) cosnd, 
n 
or finally to 
jzp| 
ar = (“- COS NU Ay M, (31) 
and ['=87a"t (32) 


Equation (31) still contains the factor 4,, which has to be taken from the 
development of the known Sommerfeld solution. (This is the only point 
where this function will be explicitly needed.) From the solution of 
Keller & Blenk (1951) and Miles (1952) it is found that 

where x is the angle between the wedge bisectrix and the incident wave 
normal. With this value, (31) becomes 


nsin uz 1(2-n) 


sinnz cosnd? My > (33) 


Equation (33) expresses the relative distances of the vortex trom the vertex 
as a function of the wedge angle (”), the incident wave direction (~), the 
incident shock strength (./,), and the direction of the vortex motion #, 
measured from the (real axis) side of the wedge. ‘The angle # still being 
unknown, (33) provides only a semi-empirical formula for o and, by (32), 
for I. Equation (33) was tested by Howard & Matthews (1955), who 
found the calculated values consistently about 10°, too low, for 8=5°. 
Similarly good results were obtained by Waldron (1954). 

It may be pointed out that even at this stage, with # still undetermined, 
it is easy to find upper bounds for quantities of interest. For instance, 
according to (32) and (33), 


which has a maximum at 


n 
sin2nd = 1 — 


An upper bound for the circulation can be found with this value of #. 


a= 2, 
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4. THE KUTTA CONDITION : 

‘The condition that the velocity remains finite at the sharp edge can be 
applied also to the single-vortex model ; it determines the direction # of the 
vortex motion, as there is only one & for which the velocity at the origin 
remains finite. It will be seen, however, that the most serious shortcomings 
of the single-vortex approximation present themselves in connection with 
the Kutta condition. ‘The previous calculations herein have been pushed 
as far as possible without making use of this requirement. 

Near the origin, the potential (18) can be developed in the following 
power series : 


—1 


The velocity remains finite for ¢=0 if the coefficient of ¢” vanishes, or, 
putting ¢,.=cexp (a) again, if 
M,A,= sinn?. (35) 
7a*to" 
Equations (29) and (35) immediately vield 
4nsin® mt = 1. (36) 
According to this equation, %* depends on the wedge angle only. ‘Typical 
values are: 


3 
Q 1/2 
90) 2/3 56-6 
180) 30 


Howard & Matthews (1955) observed for 6=5° the value }=68 ,while 


the theoretical angle according to (36) is #=88 . For larger wedge angles 
the disagreement becomes even worse. The experiments made at the 


University of Michigan with 90 -wedges (plate 1) show values of 4 of about 
32 (depending slightly on the angle of incidence) instead of 56:6. A 
critical examination of the Kutta condition for the single-vortex model is 
warranted; it will lead to an explanation of the growing difficulties for in- 
creasing wedge angle. 

For this purpose, it is even necessary to investigate the original problem 
of the vortex spiral. Let s,=at¢, be the location of the vortex sheet, 
given as a function of a suitable parameter, such as the arc length A of the 
spiral (measured from the origin). Let /°, be the circulation of the spiral 
beyond the point z,; the total circulation in the whole spiral, I’, 9 will be 
designated from now on by /,. The line density of vorticity is, with the 
present sign conventions, —d/’, dA. By superposition, the potential of the 
whole spiral flow may be expressed with the help of the potential (18) as 
‘elementary solution’ 


b=au, tA, tog = = dT. (37) 


| 
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‘The Kutta condition now can be expressed as follows : 


It is seen that (35) approximates (38) in a rough way, as parts near the origin 
(c, -- 0) can give a significant contribution to the integral in (38). 

Along the spiral, the tangential velocity is discontinuous while the 
normal velocity is continuous, that is, the real part of the potential (37) will 
have a jump while the imaginary part (the stream function) is the same on 
both sides. ‘This behaviour is assured by setting up the solution in the 
form (37) (with real J’,). If the complex potentials on the two sides of the 
spiral are called 4, and 4,, their difference at any point of the spiral is real 
and given by 

(39) 


Finally, the kinematic and dynamic conditions along the vortex sheet 
must be formulated. The condition is especially simple at the origin of 
the sheet, where the velocity normal to the spiral vanishes, and the complex 
velocities on the two sides of the sheet, w, and zw,, are parallel. In view of 
(39), equality of the pressure on both sides of the sheet can be expressed by 
the equation 


Ce. ty). (40) 
dt 


It the total circulation changes with time, it is obviously impossible to have 
both w,=0 and w,=0. 

Let the angles between the sheet and the adjacent side of the wedge at 
the origin be y, and y, on the two sides of the spiral; y,;+y.=27—/. The 
flow velocity will be zero, finite or infinite at the origin according as the 
angle y is less than, equal to, or larger than z, respectively. (It can be shown 
that the sheet will have infinite curvature at ¢, =; nevertheless, the simple 
statement above remains true.) It is now evident that one of the corners 
must form the angle ~, and the other the angle 7— 8. For, if both y, and y., 
were less than 7, both w, and w, would vanish, which is impossible for 
dl, dt+0. If one corner had an angle larger than z, the velocity would be 
infinite, and the Kutta condition thereby violated. 

Now, two cases have to be distinguished. If 8~-0, the velocity is finite 
on one side, and zero on the other side, so that, in view of (40), 


Jw,|=.\ dt (say). (41) 
On the other hand, if 8 =0 and y, =y.=7, it is impossible to have w, =0 and 
still have a simple spiral flow. 'V, a velocity together with y,=7 is 


9 


only possible if the potential behaves on one side as z" with n=2, 3,... 
These stagnation points imply a profound change in flow configurations, 
making the simple spiral impossible. | Therefore, (41) must hold for 8 
and cannot be true for 8=0. This isa first indication of an essential differ- 
ence between the spirals for zero wedge angle and for finite values of £. 


rede 
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‘This point can be elaborated further as follows. ‘he flow potential of 
the spiral, as given in(37), can be written, with the help of the Kutta condition 
(38), in the form 


+ log dT. (42) 


In this expression the integrand represents the complex potential of a flow 
element in which the vortex flow is‘ paired’ with a part of the ¢”-flow in 
such a way that the Kutta condition ts always fulfilled. “Uhe velocity at the 
origin will remain finite even if only a part of the integral in (42) is considered. 

Let the domain of integration in (42) be subdivided in two parts: a short 
length A of the spiral, beginning at the origin, and the remaining domain 
A—A. Differentiation of the latter integral will give the velocity of the 
sheet near ¢=0 with no contribution from the local sheet element; this 
is regular and behaves as C7" '. For 2=0, n= 5, this velocity 1s (in general) 
not zero; in particular it does not vanish for the spiral, as only elements 
with real ¢€, have vanishing contributions, and the sign of the velocity 
depends on the sign of the imaginary part of ¢,, which 1s always the same 
for the spiral; also, d/”, never changes its sign. On the other hand, if 8-4, 
the contribution from the integral over \— A to the velocity near the origin 
always vanishes for ¢=0, and the ‘strength’ of the zero increases with &. 

Now consider the part contributed by the integral over A, which will 
introduce the velocity jump across the sheet at ¢,=0. If 9=0, this jump 
is superimposed upon an already existing mean velocity. (It can be shown 
that for n= | the contribution from A — \ tends towards the mean velocity and 
the part from \ gives the velocity ditference in (40), for the limit A - 0.) 
In contrast, for n_. }, the contribution from the small length \ not only has 
to produce the velocity ditference but the mean velocity as well; it has to 
change the character of the How profoundly—-from a velocity of the type 
toa velocity behaving as on one side and on the other side 
of the sheet. It can be seriously doubted whether such a type of solution 
exists, 1.e., Whether a spiral solution of the form (42) can be found for n— 3. 
No such difficulties arise for 8=0, and it can be shown that the required 
type of flow around the origin is possible for a sheet which is of the form 


= const.€}", for ¢,( +im,) (The reason is that a_ velocity 

field of the tvpe c,+c,C'24+. .. can be derived from the integral over 

A— A and the spiral has to be a streamline of this meld for ¢, —- 0.) 


‘he somewhat academic question of the existence will not be discussed 
anv further, and the results of the previous considerations can be summed 
up as follows. If 9=0, the spiral induces at the origin a certain mean 
velocity which provides a speedy transport of newly created vorticity, 
stretching the sheet near the origin and reducing the density of the sheet 
such that its influence near the origin is weak. If 8 —0, on the other hand, 
the vortex elements nearest to the origin have to induce their own mean 
‘transport’ velocity, which is physically possible only by a strong accumu- 
lation of vorticity near the edge, increasing in strength with the wedge angle. 


| 
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Returning to the original question, it becomes clear why the fulfillment 
of the Kutta condition becomes a poor approximation for blunt wedges. 
Mathematically, the question is, how well does (35) approximate (38) ? 
‘The approximation is poor if a high vortex density near the origin exists 
and this must be the case for blunt wedges. 

The doubts about the existence of the idealized solution (42) for n> } 
suggest a re-examination of the phenomenon in real fluids. This will be 
based on the shadowgrams taken at the University of Michigan, where 
experiments were carried out with a pressure ratio of 1-44 of the incident 
shock. ‘hese pulses may be considered weak enough for the point of view 
adopted in this paper, namely, that the vortex generation problem is essen- 
tially an incompressible flow phenomenon. In contrast, experiments by 
Waldron (1954) with pressure ratios of 1-93 and above show supersonic flow 
near the origin combined with much more complicated vortex-generation 
patterns, which will not be discussed here. 

The shadowgram shown in plate 1 reveals a startling discrepancy 
between the classical theoretical assumption—that vorticity is generated at 
the sharp edge only—and the real fluid-behaviour. ‘The shadowgram shows 
a second vortex, smaller, but very near the edge, rotating in a direction oppo- 
site to the sense of te main vortex. A qualitative explanation of this 
‘secondary vortex separation’ is easily found. ‘The flow pattern of the 
primary vortex spiral, originating at the sharp edge, must have found along 
the side of the wedge near the vortex a deceleration, or diminishing velocity 
towards the vertex of the wedge. ‘This is connected with a pressure rise, 
and therefore the boundary layer along the wall separates and forms the 
secondary vortex. ‘Thus, secondary vortex separation is a ‘viscous’ 
effect, and the basic assumption that the flow can be treated as inviscid also 
needs re-examination. 

A series of shadowgrams of the same diffraction case as that in plate 1, 
taken at different times, show however that the flow pattern with secondary 
vortex separation is still essentially quasi-steady. Boundary-layer separ- 
ation patterns are influenced by the pressure distribution and by the Reynolds 
number. It is known that domains of Reynolds numbers exist in which 
the separation phenomenon is determined by the pressure distribution only. 
In such a domain, explicit use of the viscosity for the vortex generation 
problem is not needed, just as in the case of a sharp edge. (The dissipative 
effect of vorticity has been already neglected for the primary vortex separ- 
ation.) ‘Thus the existence of quasi-steady flow patterns with secondary 
vortex separation is understandable, but it must be borne in mind that 
these patterns might be different in a different domain of Reynolds numbers. 
The pertinent Reynolds number can be taken as J), v (v=kinematic vis- 
cosity); using for [, the estimate expressed by (32) and (31), this number 
becomes (omitting factors depending on flow geometry only) 


Re=v M22 (43) 


so that, given sufficient time (and size), all values of Re may be attained, 
with the time unit fixed by the gas properties and the shock strength. 
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‘The first two factors represent a number of the order of the molecular colli- 
sion frequency in the gas, i.e. a number of the order of 10! sec’? for atmo- 
spheric air, so that for the case represented by plate 1, Re= 10°. 

What are the consequences of the secondary vortex separation effect on 
the previous results? If the flow pattern is very nearly quasi-steady, as is 
proved by experience at least for a limited Re-domain, the dimensional 
analysis of §2 remains valid for weak shocks. The numerical factors esti- 
mated in §3 become less accurate, but no change in the orders of magnitude 
can be expected. ‘The application of the Kutta condition to the single- 
vortex model becomes obviously senseless, as the secondary vortex is very 
near to the edge and has an appreciable effect there. 

Secondary vortices become increasingly important with growing wedge 
angles. ‘hey are hardly detectable in the pictures taken by Howard & 
Matthews for 8=5. A fine picture of the phenomenon for sharp wedges 
can be found in Prandtl’s article in Handbuch der Experimental-Phystk (p.20), 
also reproduced in Goldstein’s Modern Developments in Fluid Dynamics 
(p. 40). ‘The small secondary vortex would be easily overlooked by an 
observer who is not determined to find it. 

Secondary vortices will also be strongly influenced by the motion of the 
wedge, as discussed in $2. ‘Three cases are sketched in figure 1 for the 
case 8=(); the arrow indicates the velocity 7, relative to the edge. The 
strongest secondary vortices will be found in case c. In this case, v, also 
reduces the ‘mean’ velocity at the edge which transports newly created 
vorticity away from the edge. It may be conjectured that if an inviscid 
spiral solution for 8=0 exists mathematically, it can exist only beyond a 
certain limiting value of 7, 


Vo Ve =0 Yo 


Figure 1. Vortex separation with fluid motion parallel to a sharp edge. 


‘The case represented in figure Ie is the one of particular interest in the 
theory of delta wings. It is seen that the application of the Kutta con- 
dition makes no sense for a single-vortex model, from which at most a semi- 
empirical result can be expected. 

Experiments by Michael (1955) on leading-edge separation from delta 
wings show the existence of secondary vortex separation, and even more 
complicated phenomena. For small aspect ratios, the * primary’ vortex 
appears to be broken up into two regions, one very near the tip, and a main 
vortex, so that three distinct vortices can be found. If the aspect ratio is 


increased (which corresponds to an increase of 7, in the diffraction case), 


. 
(a) (b) (c) | 
| 
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the vortices are broken up further, and numerous ‘ vortical regions’ appear. 
In these cases, agreement with the results calculated from the single-vortex 
approximation is poor. 

Finally, attention should be drawn to the uneven vortex density at the 
beginning of the sheet in plate 1; other pictures also show ‘wiggles’ and 
‘knots’ along the spiral. A detailed investigation of the vortex generation 
in a viscous fluid might show that it is a truly unsteady (not quasi-steady) 
phenomenon—possibly periodic with a frequency connected with the 
stability properties of the sheet. Fortunately, this effect hardly influences 
the overall quasi-steady character of the flow, at least in the Reynolds 
number domains observed so far. 


5. ‘THE STRUCTURE OF THE CORE OF THE SPIRAL VORTEX SHEET 
IN INCOMPRESSIBLE FLOW 

Experience indicates, and the theory presented by Howard & Matthews 
(1955) confirms, the important role of compressibility in the development 
of the core of the spiral. Nevertheless, a few remarks will be made on the 
spiral sheet in an incompressible fluid, as this question leads to a challenging 
mathematical problem of considerable theoretical interest. Previous 
treatments are due to Prandtl (1922), who gave the first formulation of the 
problem, together with some special solutions, and Kaden (1931), whose 
work represents the major basic contribution to our knowledge in this 
field. 

The formulation of the problem as given in $4 is not yet complete. 
‘The condition (40) to be fulfilled along the spiral holds in this form only 
around z,=0. For any point of the spiral, the following two conditions 
(‘ kinematic’ and ‘ dynamic’) have to be fulfilled: 

(i) the normal velocity of the spiral sheet z, must equal the normal 
velocity of the fluid; 


(ii) the pressure difference between the two sides of the sheet must vanish. 

When the pressure along the spiral is computed, the time derivative of 
the potential (space-fixed) must be expressed by the time derivative along 
the spiral, i.e., along a point moving in accordance with the kinematic con- 
dition (i). ‘The two conditions are easily formulated, but there exists a 
way of writing down both conditions directly in complex form, by use of 
Kirchhoft’s law of the vortex motion, stating that every vortex element moves 
with a velocity equal to the regular part of the total velocity at the vortex 
point. ‘The notions of this law have to be redefined for a continuous sheet. 

It is well known that the regular part of the velocity is represented, in 
case of a sheet along a smooth curve, by the mean value of the velocities on 
both sides. ‘Ihe complex ‘regular part’ of the velocity will be 


‘The meaning of a ‘vortex point’ for a continuous sheet z, has to be 
defined as a point of the sheet for which J", = constant, as an observer moving 
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with such a point will never be ‘passed’ by a vortex element. (The 
definition d/°, d\ = constant is obviously wrong, as the vortex elements can be 
stretched in the course of the sheet motion.) ‘Thus, Kirchhoft’s law, ex- 
pressed by (21) for a single vortex, can be written for a continuous sheet in 


= dz, 
-( dt = const. 


giving in a nutshell both kinematic and dynamic conditions. 
Equation (45) will be applied to ‘generalized’ quasi-steady (1.e., 
similar) flows with time dependence expressed by (16). In this case, 


the following form: 


tm 
or rts (+7) 
so that (45) can be written as 
dG 


‘Lhe left-hand side has to be derived from the complex potential (+2), putting 
dl’, =" 'dG, by tinding the complex velocities on both sides of the spiral 
and taking the mean value. ‘The computation of this quantity is naturally 
the difficult part of the problem, and (48) becomes a formidable non-linear 
singular integro-diferential equation for ¢, as a function of G. 

Kaden’s results for the spiral core can be deduced and generalized by 
using a simple approximate expression for z,, in (48), namely, 


iT, 
(+9) 
where <, (and ¢,) is now measured from the centre of the spiral. (It may 


be noted that (48) remains the same for any choice of ‘similar’ point as the 
origin.) Equation (+9) expresses the assumption that z,, is the same as the 
velocity induced by a concentrated vortex of the strength 7’, (1.e., the cir- 
culation up to s,) situated at the spiral centre. It is easily seen that this 
assumption is reasonably good if the change of J’, per.turn of the spiral is 

smail compared with 
Introduction of (49) into (48), by use of the variables defined in (46) and 
(+7), leads to the equation 
1G 


de 
A 


€,—(2m—1)GC,—. (50) 
dG 
lor the solution, put 
6, =s(G)expn(G). 
Introduction of this expression into (30) leads, upon splitting into real and 
imaginary parts, to the following two equations: 
ms? (2m —1)Gss’ =0, 


—(2m—1)s*0’ = — 


— 
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After an easy integration, the final result is 
where B is a constant of integration. ‘The connection between the absolute 
value and the argument of €, is the following: 


s= BY 2(270)-™, (52) 


For m=1, the equation of the spiral in incompressible flow would be, in 
Kaden’s approximation, s#= constant, for @—- x. 
Expressed by the radius s, it is found that 


G= (5) (53) 


1 ~ B (54) 


and on introducing a spiral are length dA for the almost circular core by the 
approximate equation 
sd, (55) 


it can be shown with the help of (52) that 
27(2m — 1)s. (56) 


The change of G per turn becomes proportional to s?, which for small values 
of s will always be less than G ~ s@”"—)™, 

The singular behaviour in the case of m= } has been pointed out before : 
the ‘concentrated vortex’ solution is exact, and there is no spiral. ‘The 
value m=} of the exponent must represent a limiting case for possible 
‘similar’ flows of this type, as, for m<}, G becomes infinite at the centre 
:=0 and decreases with s, 1.e., vorticity flows back from the core to the 
generating edge ! 


It has been noted before that for 9=0 the leading portion of the spiral 
can be handled; now approximations for the central core in incompressible 
flow have been given. The remaining task of ‘joining’ the two solutions, 
and ultimately of finding the potential for the whole flow field, will certainly 
not be easy. 
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